Density fields and halo mass functions in the Geometrical Adhesion toy Model 
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In dimension 2 and above, the Burgers dynamics, the so-called "adhesion model" in cosmology, 
can actually give rise to several dynamics in the inviscid limit. We investigate here the statistical 
properties of the density field when it is defined by a "geometrical model" associated with this 
Burgers velocity field and where the matter distribution is fully determined, at each time step, by 
geometrical constructions. Our investigations are based on a set of numerical experiments that make 
use of an improved algorithm, for which the geometrical constructions are efficient and robust. 

In this work we focus on Gaussian initial conditions with power-law power spectra of slope n in 
the range —3 < n < 1, where a self-similar evolution develops, and we compute the behavior of 
power spectra, density probability distributions and mass functions. As expected for such dynamics, 
the density power spectra show universal high-fe tails that are governed by the formation of pointlike 
masses. The two other statistical indicators however show the same qualitative properties as those 
observed for 3D gravitational clustering. In particular, the mass functions obey a Press-Schechter 
like scaling up to a very good accuracy in ID, and to a lesser extent in 2D. 

Our results suggest that the "geometrical adhesion model", whose solution is fully known at all 
times, provides a precious tool to understand some of the statistical constructions frequently used 
to study the development of mass halos in gravitational clustering. 

PACS numbers: 98.80.-k, 98.80.Bp, 98.65.-r, 47.27.Gs 



I. INTRODUCTION 

The Burgers equation [H-ll], which describes the evo- 
lution of a compressible pressureless fluid, with a non- 
zero viscosity, was first introduced as a simplified model 
of fluid turbulence, as it shares the same hydrodynami- 
cal (advective) nonlinearity and several conservation laws 
with the Navier-Stokes equation. It also appears in many 
physical problems, such as the propagation of nonlinear 
acoustic waves in non-dispersive media the study of 
disordered systems and pinned manifolds |f|, or the for- 
mation of large-scale structures in cosmology @, 0|, see 
Q for a recent review. In the cosmological context, where 
one considers the inviscid limit without external forcing, 
it is known as the "adhesion model" and it provides a 
good description of the large-scale filamentary structure 
of the cosmic web [8[. In this context, one is interested 
in the statistical properties of the dynamics, as described 
by the density and velocity fields, starting with a ran- 
dom Gaussian initial velocity and a uniform density. 
These initial conditions are expected for generic models 
of inflation of quantum fluctuations generated in the pri- 
mordial Universe and agree with the small Gaussian fluc- 
tuations observed on the cosmic microwave background. 
In the hydrodynamical context, this setup corresponds 
to "decaying Burgers turbulence" @. 

This problem has led to many studies, especially in one 
dimension (ID) and for power-law initial energy s pec tra 
(fractional Brownian motion) Eo(k) oc k n of index [77j n. 
The two ID peculiar cases of white-noise initial velocity 
(n = 0) [H, S 10l4l2l a nd Brownian motion initial velocity 
(n = -2) [lfl[ll|U5| have received much attention. For 
more general n, it is not possible to obtain full explicit 



solutions, but several properties of the dynamics are al- 
ready known [1, H, HH . In particular, for —3 < n < 1, the 
system shows a self-similar evolution as shocks merge to 
form increasingly massive objects separated by a typical 
length, L(t) - the integral scale of turbulence - that grows 
as L(t) ~ i 2 /(™+ 3 ) j while the shock mass function scales 

as ln[n(> m)] m"+ 3 at large masses @, [13, El El ■ 

In spite of these common scalings, the range — 3 < n < 1 
can be further split into two classes, as shocks are dense 
for -3 < n < — 1 but isolated for -1 < n < 1 (Tfjj . 

In higher dimensions, the situation gets more compli- 
cated as several prescriptions for the matter distribu- 
tion (again in the inviscid limit) can be associated to 
the same velocity field, governed by the Burgers equa- 
tion. They coincide over regular regions (i.e. outside of 
shocks) but they can show significantly different behav- 
iors on the shock manifold. For instance, if one uses the 
standard continuity equation mass clusters cannot frag- 
ment but they can leave shock nodes and travel along 
the shock manifold [l^, By contrast, if one uses a 
modified continuity equation, associated with a "geomet- 
rical model" for the matter distribution, thus introducing 
the Geometrical Adhesion Model (GAM), mass clusters 
are always located on shock nodes but they do not nec- 
essarily merge when they collide (in fact, collisions can 
redistribute matter over a possibly different number of 
outgoing clusters, while conserving the total momentum) 
[3, HH . The drawback of the prescription based on the 
standard continuity equation is that the latter has to 
be numerically integrated over time depriving the knowl- 
edge of the Hopf-Cole solution of the Burgers equation of 
much of its interest. By contrast, the GAM extends the 
geometrical structure of the Hopf-Cole solution to define 
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an associated matter distribution 0, 0, l2ll - l23j , through 
Legendre transforms and convex hull constructions, so 
that both the velocity and density fields can be derived 
at any time through geometrical constructions. This is a 
very convenient property, which allows faster numerical 
computations 0,[24| as well as greater analytical insights 

una. 

In all cases, beyond ID and for generic initial condi- 
tions one has to rely on numerical experiments to ob- 
tain quantitative results in those systems. The adhesion 
model has actually been studied in previous numerical 
works @,[H-i3, in ID, 2D and 3D, using the Hopf-Colc 
solution for the velocity field (however, it was not al- 
ways realized that one has to specifically complement 
the velocity field construction to unambiguously define 
the density field). 

The goal of this article is to revisit this problem, in 
ID and 2D, with the use of a novel algorithm for the 
construction of the convex hull that is more efficient and 
more robust for the construction of the convex hull tri- 
angulations (see Sect. [TV] and App. [O for details). We 
take advantage of these simulations to investigate quan- 
tities that had not been studied previously but that are 
of great interest in a cosmological context. Thus, in ad- 
dition to the mass function of shock nodes (i.e. mass ha- 
los), we also consider the probability distributions of the 
density contrast, within spherical and cubic cells, the low- 
order moments of the density distribution and the density 
power spectrum, for which there exist specific predictions 
for both the Geometrical Adhesion Model and the 3D 
gravitational dynamics. 

In particular, the reduced cumulants of the smoothed 
density contrast can serve as both a test of the accuracy 
of our numerical codes and a guide for comparison with 
the results obtained for 3D gravitational clustering by 
iV-body codes simulating the dynamics of a pressureless 
self-gravitating fluid. They are defined as 

p ~(6%)p-i' {) 

where Sr is the filtered density field at scale R (more 
precisely the filtered density-contrast field, with <5(x) = 
(p(x) — p)/p). They were shown in [28] to take a sim- 
ple form for a top hat filter and were initially derived 
for p — 3, 4, and then in [29| for all values of p, for the 
gravitational dynamics in the large-scale limit. In d di- 
mensions we have, 

i?^oo: Sf av ^(5 + ^-^(n + 3). (2) 

In the context of the adhesion model, because prior to 
shell crossings the matter field follows the Zel'dovich ap- 
proximation the reduced cumulant values take the form 
(the result has been given for the 3D Zel'dovich approx- 
imation in [lO, HJ and extended to the context of the 
adhesion model and to other dimensions in [Hj]), 

R^oo: S 3 ^^(d-n-2) (3) 



where n is the energy spectrum index defined in Eq. (|30|) 
below. The result ([3]) only holds for n < d — 3, since for 
larger n shell crossing keeps playing a role in the large- 
scale limit [l6[ . The behavior of those quantities at small 
scale is not fully understood. It has been argued that 
they should reach a constant value (at least for power-law 
spectra). This is the case in the so-called "hierarchical 
clustering models" and it was partially checked in numer- 
ical simulations, see [HJ where an explicit description of 
the small scale plateau is proposed. More precise motiva- 
tions from first principles have been put forward although 
it has never been proved explicitly that such a family of 
solutions actually exists, and even less that it was rel- 
evant in a cosmological context. The "stable-clustering 
ansatz" introduced in [33l . 34] was such an attempt, based 
on the approximation that once objects have formed they 
fully decouple from the dynamics and keep a constant 
mass and physical size. This can also be set in a broader 
multifractal description [33, |3g|. More generally, some 
physical constraints (such as the positivity of the density 
p) can be used to obtain some information on the mul- 
tifractal spectrum whence on the statistical properties 
of the density field [!?} (for instance, the coefficients S p 
can only grow or reach a constant in a small-scale highly 
nonlinear regime). However, there is no derivation of 
the precise form of the multifractal spectrum either from 
systematic approaches or well-controlled models. As de- 
scribed in this article, the Geometrical Adhesion Model 
offers the opportunity to check the validity of large-scale 
limits such as ©-(El, w hile showing a non-trivial but 
well-understood small-scale limit. 

Another focus of this paper is the Press-Schechter for- 
malism [HJ], which is widely used in cosmological large- 
scale studies. Simulations of the formation of large- 
scale structures in cosmology have indeed shown that for 
Gaussian initial conditions, such as those studied here, 
the mass function of halos defined by a given density 
threshold (typically p/p = 200) is reasonably well de- 
scribed by the Press-Schechter formula. This heuristic 
approach states that the fraction of matter, F(> M), 
that is enclosed within collapsed objects (infinitesimally 
thin shocks in the present adhesion model) of mass larger 
than M is given by the probability that, choosing a La- 
grangian point q at random, the mass M around this 
point has just collapsed to a point at the time of interest 
if ones assumes spherical collapse dynamics. For Gaus- 
sian initial conditions this reads as 

F PS (>M) = / ^-fpa(/) (4) 

Ju(M) v 

with (including the usual normalization factor 2) 

hs{v) = ^ve-^l\ (5) 
The value of v(M) can be written as 

»M = ^fy (6) 
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where a(M) is the rms density fluctuation at scale M 
and S c is determined by the dynamical evolution of the 
density field. It is 8 C w 1.69 for 3D and S c « 1.47 for 
2D, for the gravitational dynamics. In the context of the 
Burgers equation, prior to caustic formation particles fol- 
low the linear displacement field, which leads to S c = d 
for d dimensions. The factor 2, which we have inserted in 
Eq. ([5]) to ensure the correct normalization to unity, is not 
accounted for in this simple approach but can be using 
the random walk approach (as originally shown in [391 ]) 
that (at least partially) takes into account the cloud-in- 
cloud effects. It also implies that at late time all mat- 
ter points are comprised in a halo, which is expected 
to be the case for both the gravitational dynamics and 
the adhesion model. Although 3D realistic cosmological 
numerical experiments [40j-|42[ show deviations from the 
simple Press-Schechter model (|4])-([5]), the mass functions 
built by the gravitational dynamics are still very well de- 
scribed by the scaling but with a slightly different 
function f(u) than ([5]). It has been argued that those 
differences could be accounted for by various refinements 
(ellipsoidal collapse [Hj], colored noise [44[, etc.) that 
all could be similarly implemented in the GAM context. 
One of the aims of this paper is therefore to test the 
mere validity of this Press-Schechter scaling, at the level 
of the one-point mass function, within the Geometrical 
Adhesion Model, which is better controled and provides 
a much larger range of masses than 3D gravitational sim- 
ulations. This also allows us to check (especially in ID) 
the predictions that can be obtained for the large-mass 
tail of the halo mass function. 

This article is organized as follows. In Sect. [TT| we 
recall the Burgers equation and its Hopf-Cole solution 
for the velocity field. Then, we describe the associated 
"Geometrical Adhesion Model", which defines the mat- 
ter distribution that we study here. We also present the 
power-law Gaussian initial conditions that we focus on. 
They correspond to the initial conditions that appear 
in the cosmological context, for the formation of large- 
scale structures in the Universe, and they give rise to 
self-similar dynamics. We briefly present in Sect. ITni our 
numerical results for the ID case, where they can be 
checked with the help of the known analytical results. 
The large dynamical range also allows a precise test of 
scaling laws and asymptotic tails. We discuss in more 
details our results for the 2D case in Sect. IIVI After a 
brief description of our numerical algorithm, we study 
the shock mass functions that we obtain and the depen- 
dence of the low-mass and high-mass tails on the slope of 
the initial power spectrum. Then, we present our results 
for the density probability distributions and the density 
power spectrum. Next, we briefly discuss in Sect. IVlthc 
case of separable initial conditions in arbitrary dimen- 
sions, where exact results can be obtained. Finally, we 
conclude in Sect. IVII 



Note that the paper contains a few appendices that 
describe the algorithms, compare them with previous nu- 
merical studies, and present our detailed results for the 



ID and separable cases. 

The reader who is mostly interested in our results and 
the comparison with behaviors observed for 3D gravita- 
tional clustering, may skip the section IPTl devoted to the 
definition of the dynamics, and directly go to Sect. IIII1 



II. BURGERS DYNAMICS AND 
GEOMETRICAL MODEL 

A. Equation of motion and Hopf-Cole solution for 
the velocity field 

We consider the d-dimensional Burgers equation [l| in 
the inviscid limit (with d > 1), 



dtu + (u • V)u = i/Au, 



(U 



(7) 



for the velocity field u(x,£). As is well known, for curl- 
free initial velocity fields the nonlinear Burgers equa- 
tion ([7J can be solved through the Hopf-Cole trans- 
formation p5l [ill ], by making the change of variable 
ip(x,t) = 2ivlnH(x, t). This yields the linear heat equa- 
tion for S(x, i), which leads to the solution 



ipfat) = 2i/ln 



dq 



exp 



i>o (q) |x - q| : 



2;/ 



4ut 



(8) 

Then, in the inviscid limit v — > + , a steepest-descent 
method gives p], Q 



ip(x, t) = sup 
q 



-00 (q) 



|x- q 



- 1 



2t 



(9) 



If there is no shock, the maximum in ^ is reached at a 
unique point q(x, t), which is the Lagrangian coordinate 
of the particle that is located at the Eulerian position x 
at time t [H Q (hereafter, we note by the letter q the 
Lagrangian coordinates, i.e. the initial positions at t = 
of particles, and by the letter x the Eulerian coordinates 
at any time t > 0). Moreover, this particle has kept its 
initial velocity and we have 



u(x,i) = u [q(x,*)] = 



x - q(x,i) 



(10) 



If there are several degenerate solutions to © we have a 
shock at position x and the velocity is discontinuous (as 
seen from Eq. (fT0|) . as we move from one solution q to 
another one q+ when we go through x from one side of 
the shock surface to the other side). 

The solution © has a nice geometrical interpretation 
in terms of paraboloids 0, Thus, let us consider the 
family of upward paraboloids "P X! c(q) centered at x and 
of height c, with a curvature radius t, 



P*,c(q) = 



2t 



(11) 



Then, moving down P x ,c(q) from c = +00, where the 
paraboloid is everywhere well above the initial potential 
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■00 (q) (this is possible for the initial conditions (|25l) be- 
low, since we have 1-00(1)1 ~ q^ 1- ™)/ 2 , which grows more 
slowly than q 2 at large distances) , until it touches the sur- 
face defined by 0>o(q), the abscissa q of this first-contact 
point is the Lagrangian coordinate q(x, t) . If first-contact 
occurs simultaneously at several points there is a shock 
at the Eulerian location x. One can build in this manner 
the inverse Lagrangian map, x i-> q(x, t). 



form 
if(x, t) = sup 



xq- — +t0o(q) 



(15) 

Here we used the standard definition of the Legendre- 
Fenchel conjugate /*(s) of a function /(x), 



/*(s) = £ s [/(x)] =sup[s-x-/(x)]. 



(16) 



B. Geometrical Adhesion Model for the density 
field: Legendre conjugacy and convex hull 

To the velocity field u(x, t), defined by the Burgers 
equation (J7), we associate a density field p(x, t) gener- 
ated by this dynamics, starting from a uniform density 
Po at the initial time t = 0. The latter obeys the usual 
continuity equation outside of shocks. However, along 
shock lines, where the inviscid velocity field is discontin- 
uous, it is possible to define several prescriptions for the 
evolution of the matter distribution in dimensions higher 
than 1. In this article we use the prescription described 
in details in [2l|, where both the velocity and density 
fields are first defined for finite v and the inviscid limit 
is taken on a par, in a fashion which allows to derive the 
matter distribution from a geometrical construction in 
terms of convex hulls. In terms of the continuity equa- 
tion, this corresponds to adding a specific diffusive term 
that is proportional to v. Then, this term vanishes in 
the inviscid limit outside of shocks but it has a nontrivial 
nonzero limit along shocks (just as the diffusive term in 
Eq.© has a nontrivial inviscid limit, which prevents the 
formation of multistreaming flows in u(x, t)). 

We now describe how the matter distribution is ob- 
tained within this "Geometrical Adhesion Model". Let 
us first recall that an alternative description of the Burg- 
ers dynamics to the paraboloid interpretation |11|) is pro- 
vided by the Lagrangian potential 0, 0, l2l| . Thus, 
let us define the "linear" Lagrangian potential y>i(q, t) 
by 



2 



£0o (q), 



(12) 



so that in the linear regime the Lagrangian map, q i— > x, 
is given by 



*L(q,i) 



dq 



q + tu (q). 



(13) 



Thus, we recover the linear displacement field, x^ — q = 
tuo(q), which is valid before shocks appear, as seen in 
Eq.(fT0)) above. Next, introducing the function 



JJ(x,t) = L ^-+ty(x,t), 



(14) 



Therefore, Eulerian quantities, such as the velocity field 
u(x, t), which can be expressed in terms of the velocity 
potential ip(x.,t), whence of -ff(x, t), can be computed 
through the Legendre transform (fTB"]) . In particular, this 
yields the inverse Lagrangian map, x 4 q, q(x, t) be- 
ing the point where the maximum in Eq.® or (|15[) is 
reached. 

In ID, one can derive the direct Lagrangian map, 
q i y x, from this inverse map, x i— > q, using the fact 
that both maps are monotonically increasing as parti- 
cles cannot cross. However, in higher dimensions this 
is no longer the case and one must explicitly define the 
evolution of the matter distribution. As explained in 
[2lj . within an appropriate inviscid limit for the density 
field it is possible to identify the "Lagrangian-Eulerian" 
mapping q -s-> x with the Legendre conjugacy associated 
with Eq. (TT5)) . Thus, one obtains the direct map, qKx, 
by "inverting" Eq. (fT5j) through a second Legendre trans- 
form, 

y>(q, t) = £ q [i? (x, *)] = sup [q x ff(x, *)] . (17) 

x 

From standard properties of the Legendre transform, this 
only gives back the linear Lagrangian potential y>i(q, t) 
of Eq. p3|) if the latter is convex, and in the general case 
it gives its convex hull, 



ip = conv(^i). 



(18) 



Then, q and x are Legendre-conjugate coordinates and 
they are given by 



, , dH , . dip 
q(x,i ) = -s- , x(q,i) = — . 

ox aq 



(19) 



Thus, both maps, q(x) and x(q), derive from a convex 
potential and we can define the density field from these 
mappings by the conservation of matter 0, HH, H?} , 



which reads as 
Po 



p(x,t)dx = p dq, 



,1,t 'ix)= det (l 



(20) 



(21) 



the maximum Q can be written as the Legendre trans- 



Here we used the fact that both determinants are posi- 
tive, thanks to the convexity of H(x) and p(q). Thus, the 
"Lagrangian-Eulerian" mapping q -H> x ([T9"|) and the den- 
sity field (HU) define the "Geometrical Adhesion Model" 
that we study in this article. 
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Here we must note that it is possible to use other pre- 
scriptions for the evolution of the distribution of matter. 
For instance, one can use the standard continuity equa- 
tion. However, in this case one needs to numerically inte- 
grate the continuity equation over all previous times, so 
that one looses the advantages of the Hopf-Cole solution, 
which allows to integrate the dynamics to obtain at once 
the velocity field at any time through Eq.©. By con- 
trast, within the approach studied here, the density field 
at any time is fully determined by the Legendre trans- 
form (|T7)) . that is, by the convex hull (|T%|). Therefore, 
the nonlinear dynamics has been reduced to a one-time 
geometrical problem of convex analysis. In order to ob- 
tain the velocity and density fields at time t it is sufficient 
to compute the Legendre transforms (fT5j) and (fl7|) . the 
latter being equivalent to the direct computation of the 
convex hull (j 18[) . In particular, there is no need to com- 
pute the evolution of the system over previous times. 

We must point out that, although different prescrip- 
tions coincide in regular regions (outside of shocks), they 
can lead to very different behaviors on the shock mani- 
fold. Thus, it has been shown that if one uses the stan- 
dard continuity equation, limit trajectories are unique 
19, 20], so that trajectories that pass through a point 
at a given time coincide at all later times. Then, halos 
cannot fragment (particles which have coalesced remain 
together forever) but they can stop growing and leave 
shock nodes (while remaining on the shock manifold). 
In contrast, as described in details in [2l| (see also [4(), 
within the approach (|T9"| studied in this article halos can 
fragment in dimension 2 or higher. More precisely, in 2D 
the matter within shock nodes is redistributed through 
only two kinds of events, "(2 — > 2) flips" and "(3 — > 1) 
mergings" . In the first case, 2 halos collide and give rise 
to 2 new halos, whereas in the second case, 3 halos collide 
to form a single object. In each case there is a redistri- 
bution of matter but the total momentum is conserved. 
In 3D there are (2 -> 3), (3 -> 2) and (4 -> 1) events. 

On the other hand, as noticed above, in ID there are 
no ambiguities as the map q(x) is sufficient to build x(q), 
and all prescriptions coincide (in the inviscid limit). 

The "geometrical model" defined by the Legendre con- 
jugacy (IT5t - (fT7t leads to specific tessellations of the La- 
grangian q-space and the Eulerian x-space 0, [U |47| . 
More precisely, the Eulerian-space tessellation is fully de- 
fined by the Hopf-Cole solution and for the power-law 
initial conditions (|251) that we consider in this article one 
obtains a Voronoi-like tessellation. Eulerian cells corre- 
spond to empty regions (i.e. voids), which are associated 
to a single Lagrangian coordinate q as for all points x 
in a cell the maximum in ([9]) is reached for the same 
value of q. The boundaries of these cells correspond to 
shock lines in 2D (or shock surfaces in 3D) where the 
velocity field is discontinuous, and are reminiscent of the 
filaments and sheets observed in the 2D and 3D grav- 
itational dynamics. However, for the power-law initial 
conditions (|25p all the mass is contained within pointlike 
clusters located at the summits of these Voronoi-like dia- 



grams. Moreover, thanks to the geometrical construction 
that underlies the "geometrical model" (|T9l , within this 
approach this Voronoi-like tessellation is associated in a 
unique fashion to a dual Delaunay-like triangulation in 
Lagrangian space. Thus, each shock node is associated 
with a triangle in 2D (a tetrahedron in 3D) of this q- 
space triangulation, which gives the mass and the initial 
location of the particles that make up this mass cluster. 
Then, as time grows these tessellations evolve in a specific 
manner, so that these dual constructions remain valid at 
all times. This implies for instance in 2D that a colli- 
sion between two shock nodes can only give rise to two 
new shock nodes, and not to a single larger mass clus- 
ter, because two triangles cannot merge to build a single 
larger triangle (this requires three triangles embedded in 
a larger one, which corresponds to a three-body collision 
in Eulerian space) [1, [2l| ■ 

Here we may note that standard Voronoi tessellations 
have also been used in cosmology to study the large-scale 
structures of the Universe, as they provide a model of 
these large-scale structures which can reproduce some 
properties of the observed galaxy distribution |48l - l50| . 
The facts that the Burgers dynamics leads to g eneralized 
Voronoi cells as described above, see also j47|, and that 
this model provides a good description of grav itational 
clustering at large scales in cosmology [H, 0, HH , provide 
a further motivation for the use of Voronoi tessellations 
in this context. 

The fact that within the approach defined by the "geo- 
metrical model" the system can be integrated is obviously 
a great simplification. This allows both to gain a better 
understanding of its properties, taking advantage of this 
geometrical interpretation [21| , and to devise efficient nu- 
merical algorithms. This has already been investigated in 
previous works, such as ItJ. These nice properties are the 
main motivations for the use of the mapping (|19[) . rather 
than alternative prescriptions which keep the standard 
continuity equation even at shock locations but cannot 
be integrated in a similar fashion. Moreover, as we shall 
describe in this article (see also 0, [26j]) the density 
fields generated by this "geometrical model" show many 
properties that are similar to those observed in the large- 
scale structures built by the gravitational dynamics that 
is relevant in cosmology. 



C. Initial conditions 

Since there is no external forcing in Eq.([7]), the stochas- 
ticity arises from the random initial velocity Uo(x), which 
we take to be Gaussian and isotropic, whence (u) = by 
symmetry. Moreover, as is well known if the initial 
velocity is potential, Uo = — VV>0j it remains so forever, 
so that the velocity field is fully defined by its potential 
V>(x, i), or by its divergence #(x, t), through 

u = -V^, 6 = -V • u = Aif>. (22) 
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Normalizing Fourier transforms as 

0(x) = / dk e ik x 0(k), 



(23) 



the initial divergence 9q is taken as Gaussian, homoge- 
neous, and isotropic, so that it is fully described by its 
power spectrum Pg (fc) with 

(9 Q ) = 0, (0 o (ki)0o(k 2 )) = fo(ki +k 2 )F 9o (fc 1 ), (24) 

where we note <5d the Dirac distribution. In this article 
we focus on the power-law initial power spectra, 



D 



in+Z—d 



with 



3 < n < 1, (25) 



which defines the normalization D of the initial condi- 
tions. Thus, the initial conditions obey the scaling laws 

A>0: 0o(A -1 k) '= A d -<"+ 3 )/ 2 0„(k), (26) 
o (Ax) ^ A-("+ 3 '/2 o (x), (27) 



where " means that both sides have the same statis- 
tical properties. This means that there is no preferred 
scale in the system and for — 3 < n < 1 the Burgers dy- 
namics will generate a self-similar evolution [j| . This 
is why we only consider the range — 3 < n < 1 in this 
article. For the initial velocity and potential this yields 
for any A > 0, 

u (Ax) ^ A-(" +1 >/ 2 u (x), ^o(Ax) 't= w A^")/ 2 ^(x). 

(28) 

Since we have u(k,t) = i(k/fc 2 )6>(k, t), the initial energy 
spectrum is a power law, 

(iio(ki) • uo(k 2 )) - fe(ki + k 2 )£ (fci), (29) 



with Bq (A) = k~ 2 P ga (k) 



D 



(27r) d 



(30) 



whereas the initial velocity potential power spectrum 
reads as 

(^o(ki)^o(k 2 )> =fo(ki +k 2 )P^ (fc 1 ), (31) 



with P^, (fc) 



i-l-d 



(27r) d 



(32) 



point, such as the origin Xo = 0, with Uo(xo) = 0, and 
define the initial velocity in real space as 

u (x) = /"dk(e ikx -e ikxo )u (k), for -3<n<-l. 



(33) 

Alternatively, one may add an infrared cutoff and focus 
on much smaller scales (i.e. push this cutoff to infinity 
in final results). In the numerical simulations below we 
choose this second alternative as we always define our 
system on a finite box with periodic boundary conditions. 
In any case, the second-order velocity structure function 
S Uo does not suffer from this IR divergence and it reads 
as 



S , „ (xi,x 2 ) = (|u (xi) - u (x 2 )| 2 ) = DI n x~ 



(34) 



where x = |x 2 — xi|. Here we used Eq. (|3"0|) and the factor 
I n is given by 



In = 2(2tt) 



-d/2 



/ nl-d/2 

dkk n { ^ rT7 ^-k l - d l 2 J d/2 _ 1 {k) 



This reads as 

d=l : /„ = 
d = 2 : I n = 



\T(d/2) 
2sin(ri7r/2) 



r(-n) sin[(n + 1)tt] ' 
2 rl+1 sin(n7r/2) 
r[(l-n)/2] 2 sin[(n + l)7r] 



(35) 

(36) 
(37) 



Note that /„ is only defined for —3 < n < —1 as it 
diverges for n > — 1. 



2. "UV class' 1 



-1 < n < 1 



For — 1 < n < 1 the initial velocity is homogeneous but 
it is no longer a continuous function (if we do not add an 
ultraviolet cutoff). For instance, in the case {n = Q,d = 
1} it is a white noise. Thus, the initial one-point velocity 
variance, (|u | 2 ) = J dkE {k), shows an UV divergence. 
Then, it can be convenient to consider the initial velocity 
potential ipo, which is continuous (but not homogeneous), 
rather than the initial velocity Uo. Its initial second-order 
structure function reads as 

S 0o ( Xl ,x 2 ) = ([^ (x 1 )-^o(x 2 )] 2 ) =DI n _ 2 x~ n+1 , 

(38) 

where the coefficient I„_ 2 is again given by Eqs. (|3"5"|) - (l3"71) . 



1. "IR class": -3 < n < -1 

For —3 < n < —1 the initial velocity field is a contin- 
uous function but it is not homogeneous and only shows 
homogeneous increments (if we do not add an infrared 
cutoff). For instance, in the case {n = — 2, d = 1} it is 
a Brownian motion. Then, one may choose a reference 



D. Density contrast and linear mode 

In order to follow the evolution of the matter distribu- 
tion we define the density contrast, <5(x, t), by 

5(x,t) = ^. (39) 
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Then, if we linearize the equation of motion (JT)) and the 
continuity equation (which holds before the formation of 
shocks) we obtain in the inviscid limit, v — > + , 



L (k,t) = 9 (k), 5 L (k,i)=i0 o (k), 



(40) 



where the subscript L stands for the "linear" mode. 
Then, when we study the system at a finite time t > 0, 
we can as well define the initial conditions by the linear 
density field <5l(x, t), which is Gaussian, homogeneous, 
and isotropic, with a power spectrum 



P SL {k,t) = t 2 Pg (k) 



D 



(2ir) d 

and an equal-time two-point correlation 

C , 5 i (xi,X 2 ) = (£z(xi,*)(Jz,(x2,t)) 

= (2tt)^ 2 



t 2 k n+3 - 



(41) 



(fcx) d / 2 



(42) 



where a; = |x 2 — xi|. Note that for any n > — 3 the 
initial density field is homogeneous, even though the ini- 
tial velocity only shows homogeneous increments when 
-3 < n < -1. 

Since we shall study the statistical properties of the 
density field smoothed over arbitrary scales x, it is 
convenient to introduce the linear density contrast 5lt 
smoothed over spherical cells of radius r, 

&Lr= [ $&(x) = / dk~5 L (k)W(kr), (43) 



V 



with 



W(kr) 



[ ^ e ikx = 2 d / 2 r(l + d/2) 
Jv * 



Jd/2(kr) 
{kr) d / 2 



(44) 



Its variance is given by 



r(d/2) Jo 



dk k d - l P SL {k)W{kr) 



(45) 

Note that a 2 is only finite over the range — 3 < n < — 1 
if d = 1 and over —3 < n < if d = 2. For higher n it 
shows a UV divergence. 



E. Self-similarity 



For the initial conditions (|25|l that we consider in this 
paper, the rescaled initial velocity potential V>o(Aq) has 
the same probability distribution as A ( - 1 ~ Tl ' ) / 2 'i/>o(q) for 
any A > 0, when we normalize by uo(0) = and V'o(O) = 
0, as seen in Eq. ([28)) . Then, the explicit solution ([9]) gives 
the scaling laws 



i>(x,t) '= i/> (t*&x,l\ , (46) 

u(x,t) L = u(t^x,l) , (47) 

q(x,t) q(t^x,l). (48) 



This means that the dynamics is self-similar: a rescaling 
of time is statistically equivalent to a rescaling of dis- 
tances, as 



A > : t->Xt, 



A ™+3 x. 



(49) 



Thus, the system displays a hierarchical evolution as in- 
creasingly larger scales turn nonlinear. More precisely, 
since in the inviscid limit there is no preferred scale for 
the power-law initial conditions (|25j) , the only character- 
istic scale at a given time t is the so-called integral scale 
of turbulence, L(t), which is generated by the Burgers 
dynamics and grows with time as in (|49[) . Hereafter we 
choose the normalization 



L(t) = (2DT) 



2\i/(n+3) 



(50) 



where the constant D was defined in Eg. (1251) . This scale 
measures the typical distance between shocks, and it 
separates the large-scale quasi-linear regime, where the 
energy spectrum and the density power spectrum keep 
their initial power-law forms, (f30|) and (l4"Tj) , from the 
small-scale nonlinear regime, which is governed by shocks 
and pointlike masses, where the density power spectrum 
reaches the universal white-noise behavior (i.e. Ps(k,t) 
has a finite limit for k^> 1/L(t)). 

This self-similar evolution only holds for n < 1, so that 
|^o(q)| grows at larger scales, see for instance Eq. (j2"5]l . 
and n > —3, so that |^o(q)| grows more slowly than 
q 2 and the solution (j9]) is well-defined This is the 
range that we consider in this paper. The persistence of 
the initial power law at low k for the energy spectrum, 
E(k,t) oc k n+1 ~ d , that holds in such cases, is also called 
the "principle of permanence of large eddies" @ . 

In order to express the scaling law (|4^|) it is convenient 
to introduce the dimensionless scaling variables 



Q 



L(ty 



L(ty 



u 



tu 

Wr 



M 



PoL(t) d ' 



(51) 

Then, equal-time statistical quantities (such as correla- 
tions or probability distributions) written in terms of 
these variables no longer depend on time and the scale 
X = 1 is the characteristic scale of the system, associated 
with the transition from the linear to nonlinear regime. 
In particular, the variance of the smoothed linear density 
contrast introduced in Eq. (|4"5"|) writes as 



1, -3 < n < -1 : <J 2 {X/2) 



X 



-n—3 



(52) 



where /„ was given in Eq. (|3"6"|) and X — 2R (i.e. X is the 
length of the ID interval and R its radius), and 



2, -3 < n < : a 2 (R) = K n R~ 



(53) 



with 



d = 2, -3 < n< : K n = 



r(-n/2)r[(n + 3)/2] 



7r8/2(l- n )r[(l 



n)/2] 2 ' 
(54) 
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On the other hand, in terms of the dimensionless scaling 
wavenumber K and power spectra P{K) defined as 

K = L(t) k, P(K) = L(t)- d P(k, t), (55) 

the linear-regime density and velocity power spectra in- 
troduced in Eqs. BTI) and (l29l) read as 

P s (K) = - 1 - K n+3 - d , E (K) = , 1 , , K n+1 - d . 
lK ' 2(27r) d ' 01 ' 2(2n) d 

(56) 

III. THE ID CASE AS A TEST BENCH 

Our numerical implementation follows the algorithms 
described in appendix [5J At any time t, the velocity 
field u(x,t) and its potential ip(x,t) are obtained from 
the Hopf-Cole solution (fT5"j). This also gives the inverse 
Lagrangian map, x i— > q, which can be directly inverted 
to obtain the direct Lagrangian map, q 4 x, because 
both mappings are monotonically increasing. Then, the 
mapping q4x fully determines the matter distribution. 
In order to compute the Legendre transforms associated 
with Eq.(|15p we introduce the algorithm devised in [5lj . 
which first builds the convex envelope ip before taking the 
Legendre transform. This allows us to obtain the veloc- 
ity and density fields with an optimal running time that 
scales as O(N), where N is the number of grid points 
used to set up the initial conditions. By contrast, pre- 
vious works @, H3| used a slower 0(N In N) algorithm. 
We also take advantage of the fact that for the ID case 
numerous exact results are known allowing precise tests 
of the convergence properties of the codes. 

A. Shock mass function and large-mass tail 




io-° io- 4 io- a io- a o.i i io 

M 



FIG. 1: The shock mass functions N(M) obtained for sev- 
eral indices n. We plot the product M x N(M), in terms of 
the dimensionless scaling variables (|51[1 , to distinguish on the 
same plot both low-mass and high-mass regimes. The small 
error bars show the statistical error (measured from the scat- 
ter between different realizations). For n — and n = —2 
the dashed lines (which can hardly be distinguished from the 
numerical results) show the exact analytical results (|B1[) and 
(|B2| . 
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In Fig. [T]we show the shock mass function obtained 
for different values of the index n. Its shape is known 
exactly for n — and n — — 2 and the numerical results 
are shown to be in exact agreement with the theoretical 
predictions (fBl~j) and ([B2|. 

In order to measure the exponent that governs 
the high-mass tail, we show in Fig. [2l the derivative 
-dln[iV(M)]/d(M n+3 ). As shown in [tEHII, the shock 
mass function obeys the high-mass asymptotic behavior 

-3<n<l, M^oo: hxN(M) A/™+ 3 . (57) 

For — 3<n<— 2 it is possible to obtain the numerical 
prefactor using a saddle-point approach [l6| , which gives 

M™+ 3 

-3<ri<-2, M^oo; \nN(M) , (58) 

In 

where /„ was defined in Eq. (|36|) . One can check that this 
agrees with the exact result (|B2|) obtained for n = — 2 by 
different methods [3, [l5j]. For n > — 2 the relevant sad- 
dle point develops shocks and makes the analysis more 



FIG. 2: The derivative -dln[JV(M)]/d(M n+3 ) at high mass, 
for several values of n. For n = —2, —2.5, and 0, the horizon- 
tal dashed lines show the asymptotic results (|58[) and (|59l) . 
For n — —2 we also show the exact derivative obtained from 
Eq. (|B2[) (curved dashed line). For n = —1.5 the dot-dashed 
line is the value obtained from Eq. H58p . which is only approx- 
imate in this case. 

complex, although for n = it is possible to obtain an- 
alytical results and to recover the standard theoretical 
prediction (|Bip [ll|, (using our normalizations), 

M 3 

n = 0, M->oo: \nN(M) — . (59) 

We can check in Fig. [2] that each numerical curve reaches 
a constant asymptote at high mass, in agreement with 
the scaling (f57|). Moreover, for n — —2.5, —2, and 0, it is 
consistent with the analytical results ([58jl and f|59|) . 

Since these results can be recovered by a saddle-point 
approach which also applies to the gravitational 



9 



case, this suggests that the large-mass tails of the halo 
mass functions can also be exactly obtained for 3D grav- 
itational clustering, as explained in [53]. However, as 
seen in Fig. [21 the rate of convergence to the asymp- 
totic regime (|5T[) may be rather slow, especially for low n 
(but note that the deviations from the asymptotic behav- 
ior (|57[) are magnified in Fig. [5] and would appear much 
smaller in Fig. [IJ. 

More precise comparisons regarding the high-mass or 
low-mass tails can be found in the appendix [B] 

B. Press-Schechter like scaling 




i o- 4 i o-3 i o- 3 o.i i io 



With these results we are in a position to test the Press 
& Schechter formalism [38[ recalled in the introduction. 
For the Gaussian initial conditions (1251) we have for v(M), 

- 3 < n < -1 : v[M) = J A A/("+ 3 )/ 2 (60) 

V -*n 

which leads to, 



-3 < n < -1 : N } 



PS 



-l)/2 e -Af" 



_ (61) 
As noticed in [7( , it happens that Eq. (jBTj) actually recov- 
ers both the known high-mass and low-mass exponents 
of Eqs. ([57]) and ()B3j) . In fact, a saddle-point approach 
[l6| shows that it gives the exact high-mass asymptotic 
behavior ([58]) for —3 < n < —2. This is because in the 
inviscid limit i) particles move freely until shell-crossing 
in the Burgers dynamics, and therefore follow the linear 
displacement field, and ii) shell crossing only occurs for 
n > — 2 in the saddle-point that governs the high-mass 
tail ([57]). For n = —2, as noticed in [l5j . the Press- 
Schechter mass function ([6"T]) is actually exact as shown 
by the comparison with the exact Eq. (jB2l) {1-2 = 1). For 
n > — 2 the factor /„ in the exponential ([61]) no longer 
applies, as shocks come into play (note that I n actually 
diverges for n > —1), but the exponent is still valid as 
seen in ([57]) . 

In order to test whether the scaling with the re- 
duced variable v also works for the Geometrical Adhesion 
Model, we plot in Fig. [3] the function f(v) defined from 
the shock mass function as 



f(v) = MN{M) 



AM 
d In v 



2M 2 



N(M), (62) 



where we used Eg. (1601) . 

We can see that the curves obtained for n = — 1.5, — 2, 
and —2.5, are almost identical, which shows that the scal- 
ing (|62[) is a very good approximation over this range, 
—2.5 < n < —1.5, even though for —2 < n < — 1 the 
numerical factor in the exponential cutoff must show a 
weak dependence on n as explained above. In addition, 
since for n = — 2 the function f(y) defined by Eq. ([6"2"|) co- 
incides with the Press-Schechter model (JSJ, this implies 



FIG. 3: The shock mass function in terms of the scaling vari- 
able v of Eq.©, as defined by Eq.j62j. For n = -1.5, -2, 
and —2.5, the coefficient in Eq.© is given by Eq. (|36[) , but 
for n = (where 7o would diverge) it is replaced by Jo — ¥ 12. 
For n — —2 the exact curve f(v) happens to match the Press- 
Schechter model ([5]). 



that the latter is a very good approximation for this range 
of indices, —2.5 < n < —1.5. Since for n > — 1 the factor 
I n diverges we can no longer use Eq.©; this also shows 
that the scaling ()62[) can only be approximate. However, 
we display in Fig.[3]the curve obtained for n = by mak- 
ing the change Iq —> 12 in Eq.©, so as to recover the 
exact high-mass tail ([59]) . We can see that this gives a 
function f(v) that remains close to the Press-Schechter 
model (O, as could be expected from the fact that the 
latter agrees with both the exact low-mass slope (|B3[) and 
the high-mass cutoff, although a noticeable deviation can 
be seen around the peak at v ~ 1. Therefore, it appears 
that the scaling ([52"]) provides a reasonable description of 
the shock mass function for all n, provided one uses the 
appropriate normalization in the M t-¥ v relation ([6]) . 



C. Density field 

In appendix [B] we further show the behavior of the 
one-point probability distribution functions (PDF) of the 
smoothed density, and their low order cumulants, and 
compare them with known results whenever possible. 
Those tests are successfully passed. It is to be noted that 
those quantities show behaviors in qualitative agreement 
with what is expected for the gravitational dynamics, 
in both the low variance regime and the large variance 
regime, although not necessarily for the very same rea- 
sons. 

In particular, we can check in Fig. [T5] that at large 
scale the reduced cumulants S p , defined by Eq. ([B15l) (or 
the equivalent Eq.(TT])), agree with the analytical predic- 
tions (|B16]) - ([HT7] l (whence with Eq.© with d = 1), for 
n < — 2 (this upper boundary is due to strong shell cross- 
ing effects for n > —1, which are beyond the reach of 
perturbation theory). This large-scale behavior can be 
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analyzed in the same terms as for the 3D gravitational 
dynamics, through perturbation theory or saddle-point 
approaches, and both systems are similar in this respect. 

By contrast, at small scale the universal flat plateau ex- 
hibited by the reduced cumulants is a direct consequence 
of the formation of pointlike structures in the GAM. 
For the same reason, the density power spectrum (and 
also poly-spectra, although we do not explicitly show it 
here) precisely exhibits a universal k° tail in the high- 
k limit, see Fig. |T7l which is characteristic of the fact 
that formed objects are pointlike. This is not expected 
for the gravitational dynamics as the small-scale behav- 
ior of the matter spectrum, and of the reduced cumu- 
lants S p , depends on the matter profile within objects 
[37l . Then, the ratios S p do not seem to reach con- 
stant asymptotes at small scale [55j], even though their 
scale-dependence is very weak. Moreover, the "stable 
clustering ansatz" introduced in [33l . l34j ]. which would 
predict constant asymptotes (and fares reasonably well), 
is not based on such universal singularities but on very 
different arguments on the decoupling of collapsed halos, 
so that the high-Zc slope of the power spectrum depends 
on n. Thus, although both dynamics show partly similar 
behaviors at small scales, in this regime the correspon- 
dance is not exact and can be due to different physical 
processes. In spite of these limitations, some key statisti- 
cal quantities still show similar behaviors at small scales, 
such as the mass function and the PDF of the smoothed 
density described in appendix [B] Therefore, with some 
care the Geometrical Adhesion Model could still prove to 
be a useful tool to understand processes or to test approx- 
imation schemes encountered within the 3D gravitational 
dynamics. 



IV. TWO-DIMENSIONAL DYNAMICS 

We now consider the two-dimensional case, d = 2. 
Our numerical implementation follows the numerical al- 
gorithms described in appendix [Cl 

As in the ID case, in Eulerian space the velocity field 
u(x, t) and its potential i/>(x, t) are again obtained from 
the Hopf-Cole solution (|15|) . Since 2D Legendre trans- 
forms can be obtained from two successive ID partial 
Leg endre transforms, we again use the algorithm of Lucet 
51], and we obtain an optimal running time that is linear 
over the number of initial grid points, N tot = N 2 . 

The main difficulty that arises in 2D, and higher di- 
mensions, is that it is no longer possible to read the di- 
rect Lagrangian map, q m- x, whence the matter dis- 
tribution, from the "inverse" Lagrangian map, x 1— > q. 
As recalled in Sect. IIIB| within the "Geometrical Ad- 
hesion Model" this nontrivial "inversion" is performed 
through the second Legendre transform (TIT)) , or equiv- 
alent^ through the convex hull (fl8l) . Taking this Leg- 
endre transform on a grid, as in some previous works, 
would allow us to obtain an approximation of the matter 
distribution on such a grid. Then, using the same algo- 
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FIG. 4: The product M x N(M), where N(M) is the shock 
mass function in the 2D case, for several n. 

rithm as for the first Legendre transform ()15[) we would 
reach an optimal running time 0(N tot — N 2 ). How- 
ever, as explained in App. IC 21 this procedure artificially 
splits large voids into smaller voids and introduces spuri- 
ous matter concentrations. Therefore, in this article we 
prefer to exactly compute the convex hull (|18[) . without 
introducing any Eulerian grid at time t. Thus, once the 
initial conditions are given on a grid, we exactly solve 
the dynamics and we compute the exact Lagrangian and 
Eulerian space tessellations, which have been discussed 
in details in [2lj |. In particular, the density peaks are not 
restricted to a pre-defined Eulerian grid. 

Of course, the computation of the exact convex hull 
tp is a much more difficult problem than the computa- 
tion of the 2D Legendre transform on a grid. Indeed, 
whereas the latter could be reduced to ID problems, as 
explained above, the former is a standard problem of 3D 
computational geometry (since c/?(q) is embedded in 3D). 
As described in App. IG 21 and IG 41 we implement the 3D 
divide-and-conquer algorithm devised by Chan [56] . This 
recursive algorithm allows us to reach an optimal running 
time (9(iVtot In Atot)- (This is slower than the computa- 
tion of the 2D Legendre transform on a grid, because 
both problems do not have the same complexity, and the 
convex hull contains more information.) 

We discuss in more details our numerical algorithms in 
appendix[C] in particular we compare them with previous 
numerical studies in appendix IC 31 



A. Shock mass function 

We show in Fig. [4] the shock mass functions obtained 
in the 2D case (more precisely we plot the product 
M x N(M)). Thus, N(M)dM is again the mean number 
of shocks of mass in the range [M, M + dM) within a unit 
volume. As for the ID case shown in Fig.[TJ we can clearly 
see the power-law tails at low mass and the exponential- 
like cutoffs at high mass, especially for n > —1.5. For 
n = —2 and n — —2.5 where the low-mass tail grows 
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FIG. 5: The derivative -dln[iV(M)]/dln(M) at low mass. 
The horizontal dashed lines are the asymptotic behavior (|B3[) . 



faster than 1/M it is more difficult to distinguish the 
low-mass asymptote from the high-mass falloff. More- 
over, one can expect the convergence to the low and high 
mass asymptotic behaviors to be slower for lower n if the 
scaling in terms of the Press-Schechter variable v defined 
as in Eq.([S]) is still a good approximation. In agreement 
with previous works [7j, we can check that, as for the 
ID case, the shock mass function grows more slowly than 
1/M if -1 < n < 1 and faster than 1/M if -3 < n < -1. 
Therefore, the "UV class", — 1 < n < 1, again corre- 
sponds to isolated shocks, which are in finite number per 
unit volume, whereas the "IR class", — 3 < n < — 1, again 
corresponds to dense shocks, which are in infinite num- 
ber per unit volume. This agrees with the study of the 
associated Lagrangian and Eulerian space tessellations 
described in |21| . 

For —3 < n < — 1 it is again possible to obtain the 
high-mass cutoff of the shock mass function [16j . which 
now reads as 



-3 < n < -1, M -> oo 
lniV(M) 



7r -(n + 3)/2 M („ +3 )/2 )(63) 



where the factor K n was defined in Eq. ((i>4")) . For higher 
n the exponent (n + 3)/2 is expected to remain valid, 
but shocks should modify the numerical prefactor. Note 
that the analytical result (|63p now extends up to n — 
— 1, instead of n = —2 in ID. Indeed, in the general d- 
dimensional case, the overdense saddle point associated 
with this high-mass tail is only affected by shocks for 
n > d— 3. Unfortunately, the mass range of the numerical 
computations is too small to see the convergence to the 
asymptotic behavior (j63]), although they are consistent 
with the scaling M (n+3 ^/ 2 . 

At low mass .previous numerical works and heuris- 
tic arguments [3] suggest that the ID power-law tail 
(|B3[) remains valid, with the same exponent (n — 
l)/2. As in Fig. we plot in Fig. [5] the derivative 
— dln[iV(M)]/dln(Af). We can see that our numerical 
results are consistent with Eq. (|B3[) . although the numer- 



FIG. 6: The shock mass function in terms of the scaling vari- 
able v of Eg. ((65)) . as defined by Eg. ((68)) . The dotted line 
labeled "PS" is the Press-Schechter model ([5]). 



ical accuracy is not sufficient to provide a precise measure 
of the exponent (we actually get slightly higher values 
than (1 — n)/2, but this might be due to logarithmic 
prefactors or to the fact that the asymptotic regime is 
barely reached at these mass scales). As already noticed 
in Fig. |U the transition through the characteristic ex- 
ponent N(M) ~ Af _1 , which marks the divide between 
dense and isolated shocks, again appears to take place at 
n = — 1. 

Let us now investigate the Press-Schechter like scal- 
ing ((4)). In dimension d the spherical collapse relates the 
linear density contrast 5l of a spherical region to its non- 
linear density p through [la ] 



P 



(64) 



Therefore, in 2D complete collapse to a point is achieved 
at Sl = 2 and the variable v of Eq.© is now defined as 



u(M) = 



cr(M) ' 



(65) 



Then, the Press-Schechter model [38| still reads as 
Eqs.g])-©. This now yields 



2, -3 < n < : v 



where K n was defined in (|54l) . and 

n + 3 



3 < n < : Np S (M) 



, % -(n+3)/4, M (n+3)/4 

(66) 

7r -(n+5)/4 M (n-5)/4 



_ 2M („+3)/2 /(K?iX („ + 3)/2 ) 



(67) 



Thus, although the Press-Schechter prediction (|67|) re- 
covers the high-mass cutoff for —3 < n < —1, for the 
same reasons as in ID, it does not reproduce the low-mass 
tail shown in Fig. [SJ which was consistent with Eq. (|B3l) . 
This discrepancy was already noticed in previous numer- 
ical works Q. Nevertheless, it remains interesting to see 
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whether a scaling of the form (|62|) still provides a good 
approximation, albeit with a different function f(v) than 
the Press-Schechter model (J5j) . Thus, using Eq . (|66| . we 
now define f{v) as 



f{v) = MN(M) 



AM 
d In v 



4M 2 



N(M), (68) 



and we plot in Fig.[5]thc functions obtained for n < —0.5. 
We can see that the scaling by the variable v again pro- 
vides a reasonable description of the dependence on n of 
the shock mass function, although there remains a weak 
dependence on n. In particular, the low-mass exponent 
(n — l)/2 of Eq. (IB3|) corresponds to the universal low- 
v behavior f(v) oc v 2 . Thus, even though the linear 
low-f slope of Eq.© clearly fails, as seen in Fig. O a 
single quadratic slope, oc v 2 , appears to match all mass 
functions, which was not obvious a priori. However, its 
normalization shows a weak dependence on n. At high 
mass the different curves are very close, in agreement 
with (|63[) (but prefactors are expected to depend on n), 
except for the case n = —0.5 which falls somewhat be- 
low. This is expected from the constraint in Eq.(|M|. 
as for n > — 1 shocks modify the normalization of this 
high-mass asymptote. 

Note that the range of masses shown in Fig. [5] spans 
four orders of magnitude in terms of the reduced variable 
v, whereas current cosmological simulations of 3D gravi- 
tational clustering typically cover the range 0.3 < v < 4.2 
[42|, [57}, that is only one order of magnitude. This means 
that the asymptotic low-mass and high-mass tails are not 
really probed by current 3D gravitational simulations. In 
particular, they cannot measure the exponent of the low- 
v tail. For the Geometrical Adhesion Model, if the low- 
mass power-law tails (|B3[) remain valid in higher dimen- 
sions, we actually obtain jiy) ~ v at low v (with again 
no further dependence on n). This agrees with the results 
described above in ID and 2D, as well as with the separa- 
ble case in any dimension discussed in Sect. fVl below, see 
Eq. (l78l) . Current cosmological simulations cannot dis- 
criminate between such behaviors in 3D, but it would be 
interesting to check in future works whether gravitational 
clustering also gives rise to such a strong dependence on 
dimension, in terms of the reduced variable v. As seen in 
Fig. |5] for the case of the Geometrical Adhesion Model, 
such strong violations of the low-mass slope predicted by 
the simple Press-Schechter prescription do not necessar- 
ily imply strong violations of the Press-Schechter scaling 
itself. 

This is another example of the benefits that can be 
obtained by studying dynamics such as this "Geometri- 
cal Adhesion Model" , which share many properties with 
the gravitational dynamics and show complex nonlinear 
behaviors while being simple enough to provide well con- 
troled analytical and numerical analysis. They provide 
nontrivial explicit examples that can serve as a guide, 
to understand general properties or to confirm/rule out 
simple expectations. 



B. Density distribution 

We now consider the statistical properties of the 
smoothed density field. In 2D we can study the probabil- 
ity distribution function, Pr{j]) or Px{ri), of the overden- 
sity 77 within circular cells of radius r or within squares 
of size x, 



M 
^R 2 



or 77 



Pox 



£■ < 69 » 



We show both probability distributions in Fig. [7] for cells 
of the same area, that is, irR 2 = X 2 . As expected, both 
distributions are close although we can distinguish mod- 
est quantitative deviations, especially in the low-density 
tails. We can see that we recover the qualitative fea- 
tures obtained in Fig. [14] for the ID case. For —3 < n < 
— 1, the probability distributions -Pr(?7) and Px(v) show 
both high-density and low-density exponential-like cut- 
offs, whereas for —1 < n < 1 they show a low-density 
power-law tail. Moroever, for — 1 < n < 1 shock nodes 
are again isolated and in finite number, so that there is 
an additional Dirac contribution of the form PrSd(tj) or 
P x Sd{v) due to empty cells. 

At large scales we recover for — 3 < n < — 1 the Gaus- 
sian distribution associated with the linear regime [l6| . 



-3 < n < -1, R -> 00 
\nP R {Ti) ~-2 



„(n+l)/4 _ „(n+3)/4 



2R 



n+3 



K n 



/a 2 (R) 
2 



(n+l)/4 _ v (n+3)/4 



(70) 



The asymptotic behavior ([70")) holds for any finite rj if 
— 3 < n < —2, and only above a low-density threshold 
?7_, with < r\- < 1, if —2 < n < — 1 (e.g., for n = —1 we 
have rj- = 1/4). Again, for —3 < n < — 1 where typical 
density fluctuations are of order 1 7-7 — 1 1 ~ a cx _R~(™+ 3 )/ 2 , 
we can expand the argument around rj = 1 to recover the 
linear-regime Gaussian 



-3 < n < -1, R : -!> 00 

P R (V) ^ e 



\ri-l\ «i?-(»+ 3 )/ 3 : 

-(r,-l) 2 /(2<r 2 (fl))^ f 71 \ 



For — 1 < n < 0, where the linear variance (|4"5|) is still 
finite, we expect to recover the Gaussian (fTTj) at large 
scales, but the asymptotic behavior (|7D|) no longer ap- 
plies, since shocks modify the dependence on 77, whence 
the normalization of the cutoff as a function of R for any 
finite rj. For < n < 1, as for the cases — 1 < n < 1 in 
ID, where the linear variance (|45|) diverges, shocks play 
a key role at all scales and times and the density proba- 
bility distributions are always strongly non-Gaussian. 

In agreement with Fig. [71 we can expect all these 
features to remain valid for the probability distribution 
Px(v) within square cells, including the exponents in 
([70)l . but the numerical prefactors in ((70)) are modified. 

At small scales, the density probability distributions 
are again governed by the shock mass function, and the 
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FIG. 7: The probability distribution functions, Pr{vi) (solid lines) and Pxiv) (dashed lines), of the overdensity r\ within 
spherical cells of radius R and within squares of size X. The radius R is such that tvR 2 = X 2 (same cell area), with X given in 
each panel. For —1 < n < 1 there is an additional Dirac contribution (oc <5d(?7)), associated with empty cells, which does not 
appear in the figures. 



ID scaling (|B12[) becomes 



3 < n < 1, 



R 
X 



: 
■ 



Pflfa) ~ (ttR 2 ) 2 N{itR 2 V ),(72) 
P x ( V )~X*N{X*r,). (73) 



This implies in particular that the two distributions 
Pr(v) an d Px{v) coincide in the small scale limit for 
equal- area cells, in agreement with Fig. [71 However, at 
the scales shown in Fig. [7] this asymptotic regime has not 
been fully reached yet, hence we do not plot the quantity 
X A N(X 2 rj) of Eg. ([73]) to avoid overcrowding the figure. 

Using a saddle-point approach, the high-density tail of 
the probability distribution Pr(ji) reads as 



-3 < n < -1, 
In P R ( V ) ~ 



T\ — > oo : 

2?? (n+3)/2 



2R n + 



T) 



(n+3)/2 



. (74) 



which also gives rise to the high-mass tail (|63|) of the 
shock mass function. For — 1 < n < 1 shocks modify the 
asymptotic behavior, but they are expected not to change 
the exponents. Unfortunately, the range of our numerical 
computations is not sufficient to check the tails (TM1) to 
better than a factor 2, although they are consistent with 
this scaling, hxP R (rf) oc i?»+3 r? («+3)/2^ 
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n= -3.5 
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FIG. 8: The density power spectrum P(K). 



C. Density power spectrum 

We show in Fig. [8] the density power spectrum. At 
low K we again recover for all — 3 < n < 1 the linear 
regime (|56l) . whereas at high K we have the universal 
flat tail associated with shock nodes. Indeed, for the 
power-law initial conditions (|25p that we consider in this 
article, all the matter is located within point-like shock 
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nodes, which form a Voronoi-like tessellation of the Eule- 
rian space, while the boundaries of these cells are mass- 
less shock lines and the interior of the cells is empty [Il[ . 
For — 1 < n < 1, these Dirac density peaks appear to be 
isolated and in finite number per unit volume, so that 
the constant high- if tail is clearly seen in Fig. El For 
— 3 < n < — 1, shocks appear to be dense and there are 
no empty cells, in agreement with the mass functions 
and the density distributions obtained in sections IIV Al 
and IIVBI This makes the matter distribution closer to 
a continuous medium, so that the constant high- if tail 
is reached more slowly and at higher wavenumbers for 
lower n. 



V. SEPARABLE CASE IN d DIMENSIONS 

In dimensions two and higher there are no complete 
analytical results for the statistical properties of the dy- 
namics. However, it happens that the Burgers dynamics, 
and the associated Geometrical Adhesion Model, exhibit 
exact factorizable solutions in any dimension, for which 
we can obtain explicit results (especially for the cases 
n = and n = —2). This can be achieved for separable 
initial velocity potentials f58j . 

M*) = $>oW (75) 

Then, this property remains true at all times, as can be 
seen at once from the Hopf-Cole solution ©, and each 
velocity component u, (x) only depends on the coordinate 
Xi along the same direction, 

d 

7A(x,i) =Xy i) (si,*) > u i (x,t)=u®( Xi ,t), (76) 
»=i 

where the potentials ip( % > (x, t) and the velocities u^' (x, t) 
are the solutions of the ID Burgers dynamics defined 
by the initial conditions ipo \ x )- Thus, the dynamics is 
fully factorized into d ID Burgers dynamics. In terms 
of the Legendre transforms, which fully determine the 
Eulerian and Lagrangian fields as described in section [TT1 
this follows from the well-known property 

d d 

/*(s) = E/*( s *) for /(*)= !>(**)• ( 7? ) 

i=l i=l 

This states that for any function /(x) defined on R d that 
is separable (i.e., can be written as the second sum above) 
its Legendre transform /*(s) is the sum of each ID Leg- 
endre transform, as can be checked from the definition 
(|16p . This means that in d dimension, if the initial ve- 
locity potential is separable the Burgers dynamics can be 
fully factorized into d ID Burgers dynamics. This exact 
factorizability is specific to the Burgers dynamics, and 
it is not shared by more complex dynamics such as the 
gravitational or Navier-Stockes dynamics. 



As described in appendix [Dl for such factorized ini- 
tial conditions we can obtain exact results for the shock 
mass function and the density probability distributions. 
In particular, we obtain for the mass function the asymp- 
totic tails 

M -> : N(M) ~ ^/? M l, Af ( ™~ 1)/2 , (78) 
[d — 1). 

and 

M->oo: hxN(M) M (n+3)/d . (79) 

Therefore, we obtain the same asymptotic behaviors as 
those associated with the isotropic 2D case studied in 
section IIV Ai but with a logarithmic prefactor at low 
mass. This extends to any dimension d for the high-mass 
tail |16J. Thus, keeping Gaussian initial conditions with 
the scaling (|53|) preserves the characteristic exponents of 
the shock mass function, even though the isotropy of the 
system has been broken. This is not surprising for the 
high-mass tail, which can be derived from a simple saddle 
point approach and as such mostly depends on the scaling 
([53")) and the fact that the initial (linear) density field is 
Gaussian [l6[ . The robustness of the low-mass power- law 
exponent is not so obvious a priori, since it has not been 
derived in such a systematic fashion. From the analysis of 
numerical computations, the scaling (|B3[) was advocated 
using simple arguments that basically assume that the 
properties of the 2D and 3D convex hulls are similar [tJ , 
that is, governed by the scaling (|53|) . However, this also 
corresponds to assuming that the separable case studied 
in this section and the isotropic case of section IIVI give 
the same low-mass exponents, which is not obvious. 

As pointed out in Sect. IIV A[ we can note that the 
asymptotic tail ([78")) actually corresponds to a low-v 
power-law tail f(v) cx v d , in terms of the reduced variable 
v, with no further dependence on n. The large-mass tail 
([79"| also corresponds to the usual falloff f{y) ~ e _I/ I 2 at 
large v. Therefore, the Press-Schechter scaling remains 
valid in any dimension, at leading order for these sep- 
arable cases, even though the low-i' exponent strongly 
depends on the dimension. On the other hand, while the 
mass function appears more sharply peaked at higher d 
as a function of i/, as seen in Fig. 01 it flattens when it is 
drawn as a function of mass, as seen in Fig. [5] below for 
n = -2. 

The same analysis can be applied to the probability 
distributions of the smoothed density held, and we again 
recover the characteristic exponents ([74]) of the isotropic 
case. 

For the index n = — 2 we can obtain explicit expres- 
sions, which simplify in 2D as 

n = -2, d = 2: N(M) — — M~ 3 / 2 K ^2VAf) ,(80) 

P X ( V ) = H e iX rr 3/2 #o (W2 + n + ~) . (81) 

We show our results in Figs.|H]for the shock mass function 
N(M) in dimensions d = 1,2 and 3 and for the index 
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FIG. 9: The product M 2 x N(M), where N(M) is the shock 
mass function in the separable case for n = —2 and dimensions 
d — 1 (solid line), d — 2 (dashed line), and d = 3 (dot-dashed 
line). 



n = — 2. In order to emphasize the low- mass power- 
law tails we plot the product M 2 N(M) in Fig. El In 
agreement with Eqs.(|78|) and (|79[) . for higher d the mass 
function shows a smoother cutoff at high mass and a 
somewhat faster growth a low mass due to logarithmic 
prefactors. This gives more weight to extreme events, as 
is usually the case for multiplicative processes (since M = 
Yli Mi the shock mass can also be seen as the outcome 
of such a multiplicative process). This also leads to a 
broadening of the peak of the product M 2 N(M), which 
gives the fraction of matter within shocks of mass M 
per logarithmic interval of mass. We can note that such 
a flattening can also be seen in the isotropic case, by 
comparing Fig.Uwith Fig.[TJ Similar results are obtained 
for the density distribution Px(t]). 

We can note that the separable solutions studied in this 
section might serve as a basis for approximation schemes 
or perturbative expansions to describe the isotropic case 
studied in Section HVl however we shall not investigate 
further this point in this article. 



VI. CONCLUSION 

In this article we have presented a numerical analysis 
of density fields and mass functions that can be generated 
by the Burgers dynamics in the inviscid limit, the "adhe- 
sion model" in cosmology, when it is supplemented by a 
geometrical construction that explicitly defines the den- 
sity field in the shock manifolds. This leads to what we 
call the "Geometrical Adhesion Model" (GAM) for the 
density field. Our analysis focused on power-law Gaus- 
sian initial conditions, which are relevant within the cos- 
mological context, and we have considered the ID and 
2D cases. 

We furthermore have taken advantage of new more ef- 
ficient algorithms, which also make use of the geometrical 
interpretation of the system, to measure mass functions 



and density distributions over a large range of masses and 
scales. Our simulations cover seven values of the slope 
n of the initial density and velocity power spectra, that 
span the range — 3 < n < 1 where a self-similar dynamics 
develops. 

In the ID case, we have checked that our numerical 
results agree with the complete analytical results that 
are known for the two cases n — — 2 and n — 0. For 
general index n we also obtain a good agreement with 
the analytical results that apply to the tails of the mass 
function and of the density probability distributions, and 
to the low-order moments of the density contrast in the 
quasi-linear regime. In particular, this confirms the valid- 
ity of rare-event tails obtained by steepest-descent meth- 
ods. Regarding the mass functions, we found that in 
ID, they could be described with a good accuracy with 
the reduced variable, v — S c /a(M), although there re- 
mains a small dependence with n. This is the basis for 
Press-Schechter like constructions commonly used in the 
cosmological context. It also happens that the Press- 
Schechter prescription per se (e.g., derived from the ID 
spherical collapse) provides a good approximation for this 
scaling function f(v) (and it actually gives the exact mass 
function for n — —2). The density probability distribu- 
tions show the expected behaviors, with a low-density tail 
for the "UV" class — 1 < n < 1 and a sharp low-density 
cutoff for the "IR" class —3 < n < — 1. Moreover, at 
small scales we have checked that the density probability 
distributions reach their asymptotic form determined by 
the shock mass function. For the density power spectrum 
we recover the universal constant high-fc tail associated 
with shocks, which corresponds to pointlike masses in the 
density field (and discontinuities in the velocity field). 

In the 2D case we have performed a similar analy- 
sis, although the smaller range of masses and scales does 
not allow to probe with a high accuracy the rare-event 
tails. Nevertheless, our results are consistent with previ- 
ous works for the low-mass tails of the shock-node mass 
functions. We find that the scaling in terms of v still cap- 
tures most of the dependence on n of the mass functions, 
but deviations from this scaling law are slightly larger 
than in ID. Moreover, the scaling function is clearly dif- 
ferent and it shows a z/ 2 -tail at low v rather than the lin- 
ear tail obtained in ID. In this regime, as noticed in pre- 
vious works, the Press-Schechter prescription is no longer 
a good approximation. The low-density tails of the den- 
sity probability distributions show the same behavior as 
in ID, with again a qualitative difference between the 
"UV" and "IR" classes. This is related to the low-mass 
exponents of the mass functions, which are the same in 
ID and 2D (in terms of M). The density power spectra 
again reach the universal constant tail at high k due to 
the formation of shock nodes, that is pointlike masses. 

Finally, we have described how the mass functions and 
density probability distributions can be obtained in any 
dimensions for separable initial conditions, where the dy- 
namics factorizes over d one-dimensional dynamics. 

As compared with the collisionless gravitational dy- 
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namics, the nonlinear behavior of this system thus ap- 
pears as a whole much simpler to analyze as many sta- 
tistical properties can be derived from the mere fact that 
structures are all pointlike objects. As we have just seen, 
this leads to the universal flat tail for P(K) at high K. It 
also leads to constant ratios S p , defined in Eq.([T]), in the 
small-scale limit. In contrast, in the gravitational case 
relevant for cosmology (or in the Navier-Stokes dynamics 
relevant for hydrodynamics) characteristic structures are 
much more complex. Dedicated numerical simulations 
show the formation of extended halos with non-trivial 
mean density profiles and some amount of substructure 
[59M62I ] . This has prevented so far the derivation of sim- 
ple universal laws for the exponents associated with the 
density power spectrum and higher-order correlations. 

Despite these differences for the physical processes that 
take place at small scales, we find that the matter distri- 
bution generated by the Burgers dynamics, through the 
Geometrical Adhesion Model studied here, shares many 
statistical properties with the one built by gravitational 
clustering in the cosmological context. Moreover, this 
remains valid at small scales for several quantities, such 
as the mass function and the probability distributions of 
the smoothed density field. We argue then that this sys- 
tem, because of the existence of an explicit geometrical 
solution that can easily be implemented, provides a good 
tool for understanding the nonlinear processes that are 
common to both systems. One example of this is to be 
found in [63j where we explored the behavior of response 
functions (propagators) within both the Eulerian and La- 
grangian frameworks. Another line of investigation which 
remains to be explored is the use the Burgers dynamics 
as the basis of new approximation schemes, for instance 
through perturbative methods, for the 3D gravitational 
dynamics itself. 



and we use the discrete Fourier transform 

N/2-1 

u Q (x)=u (x) = ^o,fe eifci > ( A2 ) 

k=-N/2+l 

where the random complex Fourier coefficients u Q g are 

independent Gaussian variables (except for tL _^ = u 
with a variance 

w> - £ < A3 > 

We simultaneously obtain the initial potential ipo(x), us- 
ing Uo(k) — —ikipo(k) from Eq. (l22[) (and taking tio.fe = 
and ip(k) = for k — 0). Of course, with this prescrip- 
tion both the initial velocity field uq(x) and potential 
ipo{x) are homogeneous, for all n (the discretization has 
introduced both UV and IR cutoffs, for n — — 1 see |78|). 



2. Computation of the ID velocity and density 
fields 

For any time t, the velocity field u{x, t) and its poten- 
tial ip(x, t) are obtained from the Hopf-Cole solution (fT5)l 
using the algorithm described in Sect IA 31 below. 

For the ID case, we do not need to use the second Leg- 
endre transform (fTTf or p^|) to compute the Lagrangian 
map x(q) since, as both functions x(q) and q(x) are mono- 
tonically increasing, x(q) can be obtained simply by span- 
ning q(x). As a consequence, our algorithm to obtain the 
velocity and density fields has an optimal runningtime 
that scales as O(N). By contrast, previous works |7l. |24| 
used a slower O(A^lniV) algorithm. 



Appendix A: Algorithms for the ID Burgers 
dynamics 

3. Computation of a ID Legendre transform by 
1. Set up of the initial conditions building a 2D convex hull 



The system is discretized on a regular grid of N points 
with a unit step, Xi — i with i = 0, .., N— 1, with periodic 
boundary conditions. As a consequence, the analysis is 
restricted to scales x such that 1 <C x -C N, and times t 
such that L(t) <C N to avoid finite-size effects. In order 
to simplify numerical computations (e.g., for Fast Fourier 
Transforms) we choose N to be a power of 2, typically 
N = 2 23 = 8388608. We can also take p = 1 so that 
each initial "particle" i (i.e. initial grid point) carries a 
unit mass. 

To implement the initial conditions (|25p we define the 
rescaled coordinates 



To compute the Legendre transform (TT5l) we use the 
algorithm devised in [5l|, which scales linearly with N, 
taking advantage of the fact that we are given iph{q) over 
an ordered grid, qj < qj+i- Thus, we first compute the 
convex envelope if of the linear Lagrangian potential tp^. 
Then, we obtain H(x) as H{x) = C x [ipL(q)] = C-x[ i p{q)], 
using the property that the tp L and its convex envelope 
ip have the same Legendre transform. Moreover, thanks 
to the periodicity of tpo(q) the Lagrangian coordinate, 
q(x — 0), of the particle that is located at the origin 
at the time t of interest obeys —N/2 < q < N/2 — 1. 
Then, since particles do not cross each other, so that 
q(x) is monotonically increasing, to construct H(x) over 
the grid Xi = 0, 1, .., N — 1 we simply need to span </?l(<z) 
over the set of points {q(0), q(Q) + 1, .., q{0) + N - 1}. 
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1 2 3 4 5 6 7 8 9 10 q 

FIG. 10: Construction of the convex envelope <p(q) from the 
linear potential iph{q) given on a regular grid. Moving to the 
right in the (q, ip) plane we update the step-n 2D convex hull 
(p( n '(q) as we add a new data point. 




1 2 3 4 5 6 



FIG. 11: Computation of the Legendre transform H(x) from 
the convex envelope <p(q). All lines of slope x in-between the 
slopes S23 and S34 of segments (23) and (34) in the (q, ip) 
plane make first contact from below with the vertex 3. Hence 
H(x) = xq3 — ips for S23 < x < S34. 



a. 2D convex hull 

We first obtain the convex hull ip through the following 
sequential procedure. Let us assume that at step (n), 
with n > 2, we have built the convex hull p^ n \q) of 
<Pl(<?) over the first n points of this set {9(0), ..,g(0) + 
n — 1}. At this stage, tp^ n \q) is made of p points with 
2 < p < n (because of the discretization both ipr, and its 
convex hull ip are defined by a finite number of points) . 
Moving to the next step (n + 1), we add the next point 
(<?(0) + n, ifL[q(0) + n]) and going backward we remove 
if necessary the points p,p — 1, ..,p' + 1 of the previous 
convex hull (pw until its last two vertices, p' — 1 and 
p', and the new point n + 1 turn counterclockwise in the 
(q,(pL) plane. This yields the new convex hull tp( n+1 \ 
Iterating from n = 2 up to N we obtain the p vertices of 
the convex hull ip. 

This algorithm is shown in Fig. [10] at step (9). We 
have already built the convex hull associated with the 8 
points {0, 1, .., 7} and we are adding the point 8. Mov- 
ing backward we see that we must remove the vertices 7 
and 6 and the new convex hull ip^ is made of the list 
{0,1,4,5,8}. 



b. ID Legendre transform 

Second, as in [5l|, spanning the vertices j = 1, ..,f>, and 
computing the slope Sjj+i associated with the segment 
of ip(q), we note that all x such that Sj-ij < x < 
Sj.j+i have the Lagrangian coordinate q(x) = qj, which 
also yields H(x) = xqj — (p(qj). Thus, by reading the 
discrete slopes Sjj+i of the piecewise affine convex hull 
<p(q) from left to right (whence Sjj+i is monotonically 
increasing since ip is convex), we gradually "fill in" the 
values H(xi) on the grid Xi = i with i — 0,..,N, in 
order of increasing i. Clearly, both steps (computing the 
convex hull ip and next the Legendre transform H) scale 
linearly with the number of points on the grid and are 



thus optimal (HTJl . We illustrate in Fig. [TT] this second 
step, computing H(x) on the grid from the vertices of 



Appendix B: Results for the ID case 

We present here the results obtained for the ID case, 
using the algorithm presented in the previous appendix 
lAl and for initial conditions as described in the text. 



1. Distribution of shocks 

We show in Fig. [T^] the resulting distribution of shocks 
in the position-mass plane that we obtain for one real- 
ization of the Gaussian initial conditions P5|) at a given 
time. Since we use the scaling variables (j5"Tj) the statisti- 
cal properties of the output do not depend on this time. 
In particular, the typical masses (at the onset of the ex- 
ponential cutoff of the mass function) and their typical 
length-scale (e.g. the nearest-neighbor distance) are of 
order unity. In agreement with previous works 0, IToj . 
we can see that for n = 0, which is representative of the 
"UV" class, — 1 < n < 1, we obtain a finite number of 
shocks per unit length with very few high and low mass 
objects. In contrast, for n — —2, which is representative 
of the "IR" class, — 3<n<— 1, we observe a prolif- 
eration of small shocks which appear to fill all of space 
(up to the resolution of the simulation). This agrees with 
theoretical results, which show that shocks are isolated 
and in finite number per unit length for n = (ill. [l2l f oi l . 
whereas they are dense in Eulerian space for n = —2 [l3l — 
[l5| . We obtain similar figures for other indices n in both 
characteristic classes, — 1 < n < 1 and —3 < n < — 1. 
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FIG. 12: The distribution of shocks in the position-mass 
plane. Each cross corresponds to a shock, observed at a given 
time for one realization of the Gaussian initial conditions (|25p 
with n — (upper panel) and n = —2 (lower panel). The po- 
sition X and the mass M are the scaling variables (|51[) . 



2. Shock mass function 

By averaging over many realizations, and over several 
output times for each realization (thanks to the self- 
similarity of the dynamics), we can measure the shock 
mass function N(M )dM, defined as the mean number of 
shocks of mass within [M,M + dM[ over a unit- length 
interval. We have shown our results for several values 
of the index n in Fig. [T] (to avoid having a huge vertical 
range we actually plot the product M x N(M)). We can 
clearly see the power-law regime at low mass and the ex- 
ponential cutoff at high mass, with a strong dependence 
on n. We can check that our numerical results agree with 
the exact analytical results that have been obtained for 
both cases n = [Ud, 

r +ioa dsi e~ SlM 
n = 0: N(M) = 2M / =± 



-Hoo 



and n = -2 [ll[Tl], 

n = -2 : N{M) 



2ni Ai(si) 2 

dsa e S2M Ai'(s 2 ) 
2ni Ai(s 2 ) 



(Bl) 



(B2) 



FIG. 13: The derivative -d ln[iV(M)]/d ln(M) at low mass. 
The horizontal dashed lines show the asymptotic behavior 
(|F33|) . For n = —2 the curved dashed line is the exact deriva- 
tive obtained from Eq. (|B2l) . 



In agreement with Fig. [121 the shock mass function grows 
more slowly than 1/M at low mass for — 1 < n < 1, 
which leads to a finite number of shocks per unit length, 
whereas it grows faster than 1/M for —3 < n < — 1, 
which leads to an infinite number of shocks per unit 
length because of a divergent number of small shocks 
(while the total mass remains unity). 

We have checked in Fig. [5] that the high-mass tail of 
the shock mass function agrees with the analytical pre- 
dictions (|57j) - (|59p . For n — — 2 we can also check that 
our numerical result agrees reasonably well with the exact 
derivative obtained from Eq. (|B2[) . For n — —1.5, using 
the normalization given in Eq. ([58p (i.e. 7_ 1.5) appears 
to provide a reasonable approximation to the high-mass 
asymptote. This means that for —2 < n < — 1.5 shocks 
have not significantly modified the quantitative profile of 
the saddle point. It is interesting to note that the rate 
of convergence to the asymptotic regime ([57)) decreases 
with n. This is also due to the fact that the exponen- 
tial cutoff is smoother for lower n, in agreement with the 
exponent ([57| and Fig. [TJ so that the rare-event limit 
associated with these asymptotic behaviors is reached at 
higher masses for lower n. Note that the deviations from 
the asymptotic behavior ([571) are magnified in Fig. [5] and 
would appear much smaller in Fig. [T] as the exponential 
falloff is already very steep over this mass range and one 
would not distinguish in this figure the subdominant ef- 
fect of power-law prefactors. 

At low mass, previous numerical simulations and 
heuristic arguments @, E3 suggest the power-law tail 



-3<n<l, M->0: N(M) ~ M 



(n-l)/2 



(B3) 



which has onl y be en p roved rigorously for the white-noise 
case n = [ill . [l2| and the Brownian case n = — 2 
HEj . As seen in Fig. [T3l where we plot the deriva- 
tive — dln[A r (M)]/dln(M), our numerical results agree 
with the scalings (IB3|) . and for n — —2 with the full re- 
sult (|B2[) . Contrary to the high- mass tail, the rate of 
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convergence to this asymptotic behavior is roughly the 
same for all n over —3 < n < 1. 



3. Density distribution 

We now turn to the statistical properties of the 
smoothed density field. More precisely, we study the 
probability distribution function, Px(j]), of the overden- 
sity 77 within an interval of length x, 



m 
p x 



M 



(B4) 



By conservation of matter we have (77) = 1. The exact 
expression of Px(v) is again explicitly known for the two 
cases n = —2 |15| . 



Px(v) 



e 2x 7? -3/2 e -x(^+i/n) ! (B5) 



and n — [lj 



n = 0: P x (r 1 )=P x 5 D (r,) + P$(r ) ) 
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7T_ _x± f +l °° d Sl d S2 e (si+S2)X/2+(si-s 2 Y/(iX) 
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(27Ti) 
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-,oo (27Ti) 3 
e sX( V -l) + (s 1 +s 2 )X/2+(s 1 -s 2 ) 2 /(4X) 

V 

Ai(si)Ai(s 2 )Ai(si - s)Ai(s 2 - s) 

x / dre Xr Ai(r + Sl )Ai(r + s 2 ). (B8) 
Jo 

In Eq. (|H6|) the Dirac term is associated with the nonzero 
probability, P x , to have an empty interval, in agreement 
with section [Bj] and Fig. [T2] The second term P x (rf) 
is the regular part associated with non-empty intervals. 
There is no Dirac term in Eq. (IB5|) since shocks are dense 
for n = —2, as recalled in section [Bj] and seen in Fig. [T2l 
so that the probability to have an empty interval is zero. 

We show in Fig. [14] the evolution of Px(v) as we g° 
from large scales or early times (top) to small scales or 
late times (bottom), that is, from the quasi- linear regime 
to the highly nonlinear regime. We can check that our 
numerical results agree with the exact results (|B5|) and 
(|B6[) obtained for n = — 2 and n = 0. At larger scales 
we recover a probability distribution that is increasingly 
peaked around the mean, (77) = 1, whereas at smaller 
scales an intermediate power-law regime develops. This is 
similar to the behavior observed in cosmology for the den- 
sity field built by the gravitational dynamics [111 [65l - l67j , 
starting from Gaussian initial conditions such as (1251) . 



For — 3 < n < —2, at large scales and hnite 77 one goes to 
a quasi-linear regime governed by a regular saddle-point 
with 



-3 < n < -2, A — > 00 
\nP x (v) ~- 



v (n+l)/2 _ v (n+3)/2 



X n+3 



(n+l)/2 _ ?? (n+3)/2 



/(2^ 2 (A/2)) 
2 



(B9) 



For n > — 2 shocks appear as soon as t > and mod- 
ify the numerical factor /„ in Eq. (|B9p but not the main 
exponents. In particular, for n — one has [l2j 



n = 0, A — > 00, |?7-1|>A~ 



A^ 
' 12 



\V 



1 



(BIO) 



Thus we recover the large-A and large-77 exponents of 
Eq. (lB9[) . but the functional form over r\ has been mod- 
ified. Thus, whereas for —3 < n < —2 the probability 
distribution Px(v) goes to a Gaussian at large scales or 
early times, and we recover the Gaussian initial condi- 
tions (i.e. the linear regime) , this is no longer the case for 
n > —1. Indeed, as seen from Eq. (IB9p for —3 < n < —2, 
in the limit A — > 00 typical density fluctuations have 
1 7; — 1| ~ a oc X~( n+3 ^ 2 , so that we can expand the 
argument over 77 around r\ = 1. This gives 

-3<n<-2, A^oo, |r;-l| < A- (ll+3)/3 : 



Px(v) 



-( V -l) 2 /(2* 2 (X/2)) 



(Bll) 



which coincides with the Gaussian associated with the 
linear density contrast 5l- For — 2 < n < — 1 shocks 
have a modest effect on the relevant saddle point [l6j 
and we expect to recover the Gaussian (|B11[) at large 
scales, but the asymptotic behavior (IB9|) is no longer 
valid: shocks modify the dependence on ?y, whence the 
value of the exponential cutoff for any finite 77. For — 1 < 
71 < 1, where the linear density variance (1451) diverges, 
clearly one cannot recover a Gaussian such as (|B11[) at 
large scales: shocks govern the dynamics at all scales 
and the probability distribution Px (77) is always strongly 
non-Gaussian, as explicitly shown by Eq. (|B10[ ) for the 
case n — 0. We can check in Fig. [H] that the features 
associated with either case n = and n = — 2 (such 
as the power-law tail/exponential cutoff at low densities, 
the nonzero/zero probability of empty cells) extend to 
the classes — 1 < n < 1 and — 3 < n < 1 respectively. 

At small scales the probability distribution Px(v) 1S 
governed by the shock mass function, since it is domi- 
nated by the probability to have a shock of mass M = rjX 
within the cell of size A, which gives [l2j 

-3<n<l, A^O: P x {rj) - X 2 N(r)X). (B12) 



The asymptotic behavior (|B12|) holds at fixed 77, and it 
does not describe the low-density exponential cutoff that 
is always present for —3 < n < — 1 (but is repelled to 
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FIG. 14: The probability distribution, Px(r)), of the overdensity rj within intervals of length X. Smaller X (from top to bottom) 
probe deeper into the nonlinear regime. For n = —2 and n = the dashed lines are the exact analytical results (|B5|) and (|B8|) . 
For —1 < n < 1 there is an additional Dirac contribution (oc 5d(j])), associated with empty cells, that does not appear in the 
figures. In the last panel (X — 0.125) the dot-dashed lines are the asymptotic behavior (|B12I) . 



f] — > as X goes to zero). One can explicitly check 
on the exact expressions obtained for N(M) and Px(v) 
in both cases n — and n — —2 that they agree with 
(|B12|) . Since the analytical expression of the mass func- 
tion N(M) is only known for these two cases, n = 
and n = —2, we plot in the lower panel of Fig. [14] the 
asymptotic quantity X 2 N(rjX) (dot-dashed lines) ob- 
tained using the mass functions measured from our nu- 
merical computations and shown in Fig. [T] We can see 
that the behavior (|B12[) can already be clearly seen at 
X ~ 0.125, especially for -1 < n < 1. For — 3 < n< — 1 
the low-density falloff has not been repelled to very low 
r\ yet so that the convergence to (|B12[) only appears for 
f] > 1. 

In order to see more clearly the high-density tail of the 
probability distribution Px (?7) we show in Fig. [15] the 
quantity — \nP x (r])/(X n+3 ri n+3 ). Indeed, in this rare- 
event limit one can still apply a saddle-point approach, 
which yields [l6[ 



—3 < n < —2, T] — > 00 

n n+3 

\nP x {ri) 



X 



n+3 



2a 2 {X/2) 



V 



n+3 
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Indeed, this is governed by the same saddle point as the 
one associated with (|B9|) . even though we now consider 



the limit of large 77 at fixed X, instead of large X at 
fixed rj. For n > —2 the saddle point develops shocks, 
which modify the numerical factor in (|B13I) but not the 
exponents. In particular, for n = this gives 



n = 0, rj — > 00 : lnPx(?y) 



12 



(B14) 



Of course, one can check that the asymptotic behaviors 
(|B13|) and (|B14|) agree with the full expressions (|B5|) and 
(IB8|) obtained for n — —2 and n = 0. We can see 
in Fig. [15] that our numerical results reach a constant 
asymptote at high density, in agreement with the gen- 
eral scaling (|B13|) . and for n < —2 and n = they are 
consistent with the theoretical values (|B13|) and (|B14[) . 
For n = —1.5, using the normalization given in Eq. (|B13p 
(i.e. /_i.5) again appears to provide a reasonable approx- 
imation to the high-density asymptote (albeit slightly 
lower). In agreement with Fig. [5] this means that for 
—2 < n < —1.5 shocks have not significantly modified 
the quantitative profile of the high-density saddle point. 
As for the high-mass tail of the mass functions displayed 
in Fig. [5] the convergence to the asymptotic behavior 
(|B13|) is slower for lower n. We can also note that the 
shape of the function — \nP x ('rj)/(X n+3 ri n+:i ) depends 
on scale, as it typically reaches its asymptote from be- 
low at large X and from above at low X. As shown 
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FIG. 15: The ratio — In Px(r])/(X n+3 r] n+3 ), which characterizes the high-density tail of the probability distribution Px{rj). 
The horizontal dashed lines are the exact asymptotic values (|B13f) and (|B14[) for n = —2.5, —2, and 0. For n = —2 the curved 
dashed line is the exact ratio obtained from Eq. (|B5[) . For n = —1.5 the dot-dashed line is the value obtained from Eq. (|B13 |l . 
which is only approximate in this case. 



by the exact ratio obtained from Eq. (|B5l) for the case 
n = —2, which agrees with our numerical computations, 
this is not a numerical artifact. Again, note that Fig. [15] 
magnifies the deviations from (|B13j) . due for instance to 
subdominant power-law prefactors, which would not be 
easily distinguished in Fig. [14] as the exponential falloffs 
are already very steep over this density range. 



4. Low-order density cumulants 

We finally test our results with the use of the low order 
cumulants defined, S p , as 



(B15) 



They are known to reach constant values at large scales, 
and those values can be exactly computed in both the 
exact dynamics (see 16811 and references therein) and for 
the Burgers equations [79|]. Thus, for the ID case, those 
parameters reach a constant asymptote at large scales 
when —3 < n < — 1, with [16| 



-3 < n < -2, X -> oo : 
S 4 -> 



S 3 -> -3(n + l), (B16) 
3(9 + 16n+7n 2 ). (B17) 



For — 2 < n < — 1 shocks modify the large-scale asymp- 
totes, while for — 1 < n < 1 the coefficients S p typically 
vanish for odd p and go to infinity for even p, as can 
be explicitly checked in the case n = [12] where ex- 
act results can be derived from Eq. (IB8[) . On the other 
hand, at small scales they reach constant asymptotes, 
for all n in the range — 3 < n < 1, as the density cu- 
mulants are governed by the pointlike masses associated 
with shocks. The exact values of these small-scale asymp- 
totes, associated with the highly nonlinear regime and 
governed by the shock mass function, are only known for 
the two cases n = [l2| and n = — 2 [lfij]. In the case 
n = —2 it happens that the coefficients S p are actually 
scale-independent, so that the quasi-linear values (IF317|) 
hold for all X and we have UM 



S p = (2p-3)! 



(B18) 



We can check in Fig. [111 where we plot the coefficients 
S3 and 54, that our numerical computations agree with 
the results recalled above. In particular, we clearly see 
the small-scale universal constant asymptotes, due to 
shocks, except for the case n = —2.5. There it is not 
clear whether the deviation is due to the finite numerical 
resolution or to the slow convergence to the small-scale 
asymptote. At large scale we can see the same curve ap- 
proach the asymptote (|B17[) . As for the high- mass and 
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FIG. 16: The low-order coefficients S3 (upper panel) and 5*4 
(lower panel) denned in Eq. (|B15 [). For n — and n = — 2 the 
dashed curves are the exact analytical results. Both 53 and 
ft are constant in the case n — —2. For n = —2.5 the dashed 
line at X > 1 shows the large-scale asymptote (|B17|) . For 
n — —1.5 the dot-dashed line at X > 1 is the value obtained 
from (|B17|) . which is only approximate in this case. 



high-density tails of the shock mass function and of the 
density distribution, for the case n — —1.5 the value given 
by Eq. (|B17j) is a very good approximation for S3 at large 
scales. This again means that shocks do not have a sig- 
nificant effect for —2 < n < —1.5. For S4 the error bars 
are too large to draw any conclusion on the accuracy of 
(1B171) . 



5. The density power spectrum 

We show in Fig. [17] the density power spectrum as a 
function of the wave number K. At low K in the "lin- 
ear regime" we recover the initial condition ([56| . Note 
that this holds for all — 3 < n < 1, even though shocks 
cannot be ignored even at a qualitative level for n > — 1 
at all scales (in particular they make the real-space vari- 
ance (?? 2 )c finite even though the linear variance a 2 was 
infinite). At high K, in the highly nonlinear regime, we 
recover the universal constant asymptote, due to shocks 
(as the haloes that form are point-like it is easy to see 
that white-noise tails are expected to develop at large 
K.) For n = and n = — 2 we also plot the analytical 





\ n = -3.5 














- 


-3 V"'* 








~ 1 / 1 It 

! ^0.5/ J 0,5 | 







10- 3 lo- 2 0.1 1 10 10 a 10-' 

K 



FIG. 17: The density power spectrum P(K). For n = and 
n = —2 the dashed lines are the exact analytical results, which 
are almost completely masked by the numerical results. 

results [HI, EH as dashed lines. They reproduce exactly 
the numerical results showing that the numerical results 
are free of finite volume effects. 



Appendix C: Algorithms for the 2D Geometrical 
Adhesion Model 

1. Initial conditions and Eulerian velocity field 

As for the ID case, we discretize the system on a reg- 
ular N x N grid, with unit steps and periodic boundary 
conditions, and we choose N to be a power of 2, typically 
N = 2 11 = 2048. To implement the Gaussian initial con- 
ditions we also introduce rescaled coordinates as in (|A1[) . 
It is convenient to use the velocity potential i/'Oj which 
is obtained from a discrete Fourier transform as in (|A2jl . 
with now 

(l< t l 2 > = (^(|)""«^ (d) 

This yields the initial velocity field uo(q) through 
Eq-dlU). 

In Eulerian space, the velocity field u(x, t) and its po- 
tential i/)(x, t) are again obtained from the Hopf-Cole 
solution (fl"5|) . Thus, we need to compute the 2D Leg- 
endre transform £ x [(/?^(q, t)], over the regular N x N 
grid Xjj t i 2 , from the periodic linear potential ip l (q, t) 
defined over the regular N x N grid q^,^- Thanks to 
the period N of the system the 2D Eulerian coordinates 
x«i,«2 = (11,12), with < i < N — 1, are associated 
with Lagrangian coordinates q, lj2 = (31,32) that obey 
-N/2 < j < N + N/2 - 2. Therefore, we first extend 
ip L (q, t) to the larger grid -N/2 < j < N + A^/2 - 2 (us- 
ing periodicity) and next we compute £ x [(/3i(q, t)], using 
the well-known property that a 2D Legendre transform 
can be obtained from two successive partial ID Legendre 
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transforms, 

£x 1 ,x 2 l l PL{qi,q2)} = max[xi9i +x 2 q2 ~ ViCffi, QSj)] 

9l,«2 



max 

91 



max [x 2 92 

92 



<PL (91,92)] 



= £ X1 [-^(Vi)] 



(C2) 



Thus, for each ID Legendre transform in Eq. (|C2[) we use 
the algorithm of Lucet [5l| , used in appendix |B] for the 
one-dimensional case and described in appendix [51 tak- 
ing advantage of the fact that all functions are given on 
regular grids. This allows to compute the 2D Legendre 
transform H(x.) in linear time over the total number of 
grid points, N to t — N 2 , which is thus optimal. 



2. Direct Lagrangian map 

In addition to the Eulerian fields u(x, t) and ip(x.,t), 
the procedure (|C2[) yields the inverse Lagrangian map 
x4q. However, contrary to the ID case this is no longer 
sufficient to obtain the direct Lagrangian map q^x. In 
fact, as recalled in section Hi B[ the latter (and the asso- 
ciated density field) depends on the precise definition of 
the inviscid limit. In this article we consider the proce- 
dure described in section III B[ where the Lagrangian map 
x(q) is defined from the convex hull <ys(q) by Eq. ([T5)) . As 
shown in [21| , this corresponds to adding a specific diffu- 
sive term to the right hand side of the equation of conti- 
nuity, which vanishes in the inviscid limit v — > + except 
along shocks. As described in [2l|, see also 0,0,113]) m 
two dimensions the convex hull y(q) defines a triangula- 
tion of the Lagrangian q-space (because the convex hull 
of <^z,(q) is made of a set of triangular facets) which is 
associated with a Voronoi-like tessellation of the Eulerian 
x-space. 

From the Legendre duality (fTT)) we can see that within 
this prescription the direct Lagrangian map, q^>x, can 
be obtained from the Legendre transform of the Eulerian 
function £f(x). Thus, the position x of the particle of 
Lagrangian coordinate q is given by the point x where 
the maximum in Eq. (|17[) is reached. Therefore, since 
we have already obtained H (x) through a 2D Legendre 
transform, as explained above, we could use the same 
algorithm to apply a second 2D Legendre transform to 
-ff(x). This would give y(q) on a regular grid, as well 
as the direct Lagrangian map, q 1— > x. As noticed in 
[5lj . this procedure, based on two successive Legendre 
transforms, provides a very fast algorithm to compute on 
a grid the convex hull y(q) of any function (^(q), since 
it scales linearly with the total number of points iVtot of 
the grid (as we have recalled above for the computation 
of#(x)). 

In contrast, it is known that the explicit computation 
of the 3D convex hull scales at least as iV to t l n -Wtot • The 
reason for this longer execution time is that by "explicit 
computation of the 3D convex hull" we mean that, given 
the initial function i pl{ ( \) on a grid of A to t points q 3 , 



which defines a set of 3D points (91, 92, ^Pl)j, we want 
to obtain the subset of N v vertices that belong to the 
lower convex hull as well as its Nf triangular facets (which 
specifies how the vertices are gathered into triplets, to 
form these facets; note that each vertex can be a summit 
of several facets) . Clearly, this involves more information 
than the mere knowledge of the values of y>(q) on a grid, 
which explains the different scalings with N to t of these 
two problems (in particular, once we know the facets of 
the convex hull it is immediate to compute y(q) on any 
grid, while the converse is not true). Note that this is a 
truly three-dimensional problem. 

In spite of the explicit expression (jf 7[) . which gives the 
direct Lagrangian map x(q) through the Legendre trans- 
form of H(x), we use in this article the explicit computa- 
tion of the 3D convex hull (i.e. we compute the list of its 
triangular facets) to obtain the direct Lagrangian map, 
q 1— > x. This is necessary to obtain the Lagrangian and 
Eulerian space tessellations and to follow the intricate 
dynamics of shock nodes, which undergo both merging 
and fragmentation events. These processes are described 
in detail in [2lj . where we used for numerical compu- 
tations the algorithm that we describe in appendix IC 41 
in the present paper. If we only require snapshots of 
the Lagrangian map and of the density field, as in the 
present article, the faster Legendre transform (fT7|) would 
be sufficient as noticed above. However, in practice it in- 
troduces an additional source of error in numerical com- 
putations. Indeed, the function £f(x) being defined as 
the Legendre transform (|I5p and (fL (q) being defined on 
a set of discrete points, it is piece wise afiine. In fact, for 
the self-similar initial conditions (1251) this is not a nu- 
merical artifact due to the discretization and the planar 
facets of -ff (x), which define the Eulerian-space Voronoi- 
like tessellation, correspond to voids (i.e. empty regions) 
in the Eulerian density field. Moreover, their typical size 
scales with time as the characteristic scale L(t) defined in 
Ea. flST))) . However, if we compute -ff(x) on a grid, using 
the Legendre transform algorithm described above, it is 
clear that because of the finite precision of real numbers 
in computers such planar facets will show small wrinkles. 
Then, when we apply a second Legendre transform to 
H(x) to compute ip(q), we artificially split large voids 
into smaller voids and introduce spurious matter con- 
centrations (associated with the contact points of these 
wrinkles with their convex envelope) j§0] . This is not nec- 
essarily a serious problem if one is not interested in the 
properties of the Lagrangian and Eulerian space tessella- 
tions, as long as one makes sure to discard these spurious 
low-mass shock nodes. However, to be fully consistent 
with our previous work [2lj and to avoid introducing un- 
necessary sources of numerical error we prefer not to use 
this simple approach and to explicitly compute the 3D 
lower convex hull ip as well as the Lagrangian and Eule- 
rian space tessellations. 

Therefore, to obtain ip (more precisely, the list of its 
triangular facets) from the values (pL (q? ) of the linear La- 
grangian potential on the 2D grid of N to t points q-, , with 
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^Vtot = N 2 , we need a 3D algorithm, which computes the 
convex hull of a finite set of points in three dimensions. A 
minor simplification is that we only require its lower part, 
since ip is the lower convex envelope of ifL- A brute- force 
method, testing each triple of points as a possible facet, 
takes a running time 0(N^ ot ), whereas a slight modifi- 
cation improves this to 0(N^ ot ) by testing each pair of 
points as a possible edge [69(. Gift-wrapping algorithms 
[70l [7l| scale as 0(N tot Nf) (where N { is the number of 
output facets) by generating facets one at a time via im- 
plicit searches. Incremental methods gradually update 
the convex hull as initial points are inserted one after 
the other and can achieve an optimal 0(N to t In iVtot) 
expected running time [72| by randomizing the inser- 
tion order (so as to beat the 0(N 2 ot ) worst-case perfor- 
mance). Finally, the divide-and-conquer method, pro- 
posed in [73| for 2D Voronoi diagrams, and in [74| for 3D 
convex hulls was the earliest algorithm to achieve opti- 
mal 0(N to t In iVtot) running time. As its name suggests, 
this algorithm divides the point set into two halves by 
a vertical plane, recursively computes the hull of each 
half, and merges the two hulls into one. As usual, it is 
the bisection by two and the recursivity that bring down 
a factor A to t to IniVtot in the running time. Since we 
need to compute many convex hulls ip, with typically 
JVtot = N 2 = 2 22 = 4 1 94 3 04, in order to accumulate a 
sufficiently large number of realizations and output times 
to measure the statistical properties of the dynamics over 
several regimes, it is necessary to use a fast algorithm 
that scales as 0(N tot lnA^t)- 

We choose the 3D divide-and-conquer algorithm as im- 
plemented by Chan [56]. This provides a transparent in- 
terpretation of the method which is well-suited to our 
case, where the initial points are on a regular grid and 
we only need the lower part of the 3D convex hull. We 
describe this recursive algorithm in appendix IC 41 This 
gives the triangular facets of ip as well as the Lagrangian- 
space triangulation at any time t. Moreover, the slope 
(x\ 1 X2) of each facet gives the Eulerian-space position 
x of the shock node which contains all the matter asso- 
ciated with this Lagrangian-space triangle, with a mass 
equal to the triangle area (setting again po = 1). Then, 
listing for instance in clockwise order the facets that have 
a common vertex q c we obtain the Voronoi-like cell as- 
sociated with q c , each of these facets giving a summit 
Xj of this Eulerian-space cell. These summits are shock 
nodes whereas the cell itself is a void (i.e. with zero mat- 
ter density) and the cell boundaries are zero-mass shock 
lines. 



3. Comparison with previous numerical studies 

To compare with previous works^ let us first note that 
some previous numerical studies d, [75|, [76| of the "ad- 
hesion model" are not based on the Legendre transforms 
of Sect. Ill Bl but on the standard continuity equation. 
Thus, keeping the viscosity v small but finite, they first 



compute the velocity field at all times through the Hopf- 
Cole solution and next integrate particle trajectories 
over time using this velocity field. As we have recalled 
in Sect. Ill Bt the dynamics obtained in this fashion is in 
fact not identical to the one studied in this article (and 
some other works), as in the inviscid limit the behaviors 
of particles on the shock manifold are different. 

Other numerical works [25T - [27j have taken advantage 
of the geometrical interpretation (fTTj) . in terms of first- 
contact paraboloids, of the Hopf-Cole solution in the in- 
viscid limit to avoid integrations over time. Thus, by 
spanning the Lagrangian q— space with paraboloids (fTTj) 
of Gaussian curvature 1/t 2 at the apex, one can sepa- 
rate Lagrangian coordinates into "stuck" and "free" par- 
ticles. "Free" particles are such that the paraboloid that 
makes contact with the surface V'o(q) a t position q has 
no other intersections with the initial velocity potential 
ipo (q). Therefore, such particles have not experienced 
any shock yet, and their Eulerian location x is given by 
the apex of this first-contact paraboloid. "Stuck" parti- 
cles are such that this paraboloid has other intersections 
with tpo, which means that they have already experienced 
shocks (and their Eulerian location is no longer given by 
the apex of this paraboloid). As recalled in Sect. Ill Bl this 
geometrical construction in terms of paraboloids, which 
only relies on the Hopf-Cole solution ©, applies to all 
prescriptions (i.e. using the standard continuity equa- 
tion as well as using the "geometrical model" (|T9"j) ). It 
is sufficient to describe regular regions (i.e. outside of 
shocks) associated with "free" particles, where there is a 
one-to-one mapping q O x. In our case, for the power- 
law initial conditions ((23)) . this gives the Voronoi-like di- 
agrams associated with empty regions and their bound- 
aries, see for instance Fig. 7 in [21] and @, [25l - |27l . |47j . 
Next, scanning the Eulerian x— space with paraboloids, 
one obtains the Eulerian location of the particles that 
form the boundaries of the "stuck" Lagrangian domains. 
Moreover, considering the paraboloids that have three si- 
multaneous contact points all the particles located within 
the associated Lagrangian triangle are set to the Eule- 
rian location given by the apex of the paraboloid. In 
this fashion one reconstructs the Lagrangian-space tri- 
angulation associated with the "geometrical model" de- 
scribed in Sect. Ill Bl without performing Legendre trans- 
forms and convex hull constructions. However, this pro- 
cedure is rather intricate and the successive scans of the 
Lagrangian and Eulerian grids by paraboloids introduce 
some numerical inaccuracies. 

Finally, the use of Legendre transforms and convex hull 
constructions was introduced in [7|, [24| , from the defini- 
tion of the "geometrical model" described in Sect. Ill Bl 
This provides a very elegant method to obtain the direct 
Lagrangian map, q 1— > x, and the associated matter dis- 
tribution, without looking for first-contact paraboloids 
and trying to invert the x H> q map. Then, 0, [24| first 
compute the function H(x) through the Legendre trans- 
form (|15p . which gives the inverse map, x 1— ¥ q, and the 
Eulerian velocity field u(x). Next, they obtain the direct 
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map, q i-> x, and the associated density field, from the 
second Eq. (|19l) . by computing the nonlinear Lagrangian 
potential tp(q) from -ff(x) through a second Legendre 
transform (|17l) . 

As explained above, our procedure differs in two re- 
spects. For the first step, we also compute i?(x) and 
the velocity field from the Legendre transform (TTBl . but 
as described in Sect. IC II we use a faster algorithm that 
scales as 0(N 2 ) : instead of 0[(N\nN) 2 } in 0,H3. For 
the second step, as described in Sect. lC 2l we do not com- 
pute tp(q) through the second Legendre transform (flTl) . 
to avoid inaccuracies associated with numerical wrinkles 
in the planar facets of H(x). Rather, we directly compute 
tp as the convex hull of the linear Lagrangian potential 
ipL, from Eq. (|18p . without reference to an intermediate 
Eulerian grid. Note that this is a standard problem of 
computational geometry, that we solve exactly using a 
convex hull algorithm. Therefore, once the initial veloc- 
ity potential ^o(q) is given on a grid, whence the lin- 
ear Lagrangian potential ipL, no further sources of error 
are introduced by the use of an Eulerian grid, since we 
compute the exact convex hull tp under the form of a 
list of triangular facets whose slopes give the exact Eule- 
rian locations of the shock nodes (without reference to an 
Eulerian grid). Thus, both the Eulerian-space Voronoi- 
like tessellation and the Lagrangian-space Delaunay-like 
triangulation are exactly computed, without introducing 
further inaccuracies (up to the precision with which real 
numbers are represented by the computer). In particu- 
lar, this removes any source of ambiguities. As noticed 
in Sect. IC 21 although using the second Legendre trans- 
form (fTT|) would provide a faster algorithm that scales as 
0(N 2 ), the exact convex hull algorithm that we use in 
this article scales as 0(N 2 \nN), which is still faster than 
the method used in previous works. 



4. Computation of the 3D lower convex hull tp 

We describe here the 3D divide-and-conquer algorithm 
that we use to compute the convex hull tp (more precisely, 
the list of its triangular facets) following the implemen- 
tation introduced by Chan 56]. The main point is to 
transform the 3D problem into a "kinetic" 2D problem, 
which is easier to visualize. This proceeds by scanning 
the convex hull tp in order of increasing slope x 2 along 
direction 2, so that one only needs to study the evolution 
with "time" x 2 of a curve z(q%), as explained below. 

The input is the set of 3D points {(qi, Q2, fL)j} with 
j = 1, .., N ta t- Let us choose a given slope x 2 along the 
direction q 2 and consider the points {(qi,z)j} in the 2D 
plane (q\, z) with z = tpi — x 2 q 2 . Thus, Zj is the signed 
vertical distance between the 3D point (J) and the plane 
P X2 of equation tpi = x 2 q 2 . Then, the convex envelope 
of the set of 2D points {{qi,z)j} gives the vertices of tp 
that come into contact with planes of slopes {x\, x 2 ) (i.e. 
of equation tp^ = xiqi + x 2 q 2 + c). As one goes from 
left to right, i.e. in order of increasing q\, the slope Xi 
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FIG. 18: The projection of the 3D points {(qi, q%, <pi,)j} onto 
the (qi,z) plane. The convex hull of these points gives the 
vertices of the 3D convex hull tp observed from planes of a 
given slope X2 along direction q^. As xi increases the points 
move in the plane (qi, z) along verticals, with a "speed" equal 
to — Q2, shown by the arrows for points A, B and C. For a 
slightly larger value of X2 the point C will move below the 
segment [AB] and appear as a new vertex in the convex enve- 
lope z c (qi). This means that points {A, B, C} are a triangular 
facet of the 3D convex hull tp. 



along the axis qi also increases (as for the ID case studied 
in section 1X1]) . This is shown in Fig. fH] The N data 
points with the same coordinate q\ on the initial N x N 
grid appear on the same vertical line in the z) plane. 
By going from x 2 — -oo up to 12 = +00 one scans 
in this fashion all vertices of the lower convex hull tp. 
Therefore, one obtains a movie, running with "time" x 2 , 
where points of fixed abscissa q\ and varying height z — 
tpL — x 2 q 2 move along verticals at different speeds —q 2 , 
so that the convex envelope in the (qi,z) plane evolves 
with time, as its vertices can appear and disappear. 

In addition, one obtains the triangular facets by 
recording the insertion and deletion events. In an in- 
sertion event, a new vertex C of abscissa appears at a 
"time" x 2 in the convex hull z c (qi), in-between vertices 
A and B (this occurs when the point C crosses from 
above the segment [AB] in the plane (qi,z)). Then, the 
triplet {A,B,C} is a planar facet of tp (with a slope x 2 
along direction 2). In a deletion event, a vertex C located 
between vertices A and B of the convex hull z c (qi) disap- 
pears, and this again means that the triplet {A, B, C} is 
a planar facet of tp. In our case, where the initial points 
are located on a regular N x N grid over the (qi,q 2 ) 
plane, we also have exchange events, as a new vertex C 
can replace an older vertex C" that has the same coordi- 
nate qi . At the crossing time x 2 these vertices are located 
between two vertices A and B and we obtain two facets, 
{A, C, C'} and {B, C, C'}, of the convex hull tp. We show 
in Fig. 1181 a configuration just before an insertion event, 
as the vertex C will soon move below the segment [AB] 
(the arrows show the "velocities" —q 2 of these points). 
By computing analytically the successive crossing times 
x 2 we move from one event to the next one. Therefore, 
we do not need any discretization over x 2 and we obtain 
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tions, which fully define the distribution of matter at a 
given time, see [21| for detailed descriptions. 



Appendix D: Separable case in d dimensions 

We describe here the factorizable solutions of the dy- 
namics presented in Sect. [V] 



1. General index n 



FIG. 19: Merging of the 2D convex hulls L and R associated 
with the left and right halves of the points in the (qi, z) plane. 
The convex envelope of all the points is obtained by determin- 
ing the bridge between vertices u and v of the left and right 
hulls. As X2 increases, changes to L and R on outer sides of 
the bridge are recorded whereas changes to L and R within 
]u, v[ do not contribute. In addition, as neighbors of u and v 
cross the line (uv) the bridge is modified. 



the exact list of the facets of the convex hull p, in order 
of increasing slope X2- This fully defines p. 

The algorithm as described so far takes a running time 
0(7V t 2 ot ). In order to achieve a faster 0(N tot In iVtot) per- 
formance, we adapt the divide-and-conquer method in- 
troduced in [56| . Since the data points are given on a 
regular grid over the ((71,92) plane, the columns along 
<72 at fixed q\ (which appear as vertical columns in the 
(gi , z) plane as in Fig. [T8|) are stored in increasing order of 
qi . Then, we recursively create a movie for the left lower 
hull L of the N/2 left columns and a movie for the right 
lower hull R of the N/2 right columns. Next, we build 
a movie for the lower hull of all columns by a merging 
process. This is done by identifying the common tangent 
uv, called the bridge, and removing the part of L to the 
right of u and the part of R to the left of v. Then, start- 
ing from X2 — —00, as "time" progresses events that take 
place on either side of the bridge are recorded, whereas 
events that take place in-between vertices u and v do 
not contribute. In addition, by paying attention to the 
neighbors of u and v we update the position of the bridge. 

As for the 2D Legendre transform (fT5|) used to com- 
pute £f(x) and the velocity field, we first extend </>i(q) 
to the larger grid -N/2 <j<N + N/2 - 2 (using peri- 
odicity) and we compute the convex hull associated with 
this set of points. This ensures that boundary effects are 
eliminated for the points in the range < j < N — 1 
we are interested in. We obtain in this fashion the tri- 
angular facets of tp over this domain, which defines the 
Lagrangian-space triangulation at a given time t. More- 
over, the slope (x\,X2) of each facet gives the Eulerian- 
space position x of the shock node which contains all the 
matter associated with this Lagrangian-space triangle, 
with a mass equal to po times the triangle area. Thus, 
from the list of the triangular facets of p we obtain at 
once both the Lagrangian and Eulerian space tessella- 



For the factorizable initial conditions presented in 
Sect. El defined by Eg. ([75)1 with independent ID Gaus- 
sian fields along the different directions, the mass M 
of a shock node is the product of the "ID masses" Mj 
along directions i (since all directions are described by 
the same index n we can work with the dimensionless 
scaling variables as defined in (|5ip , with a unique char- 
acteristic length L(t) given by Eq. ([50]) ). Then, the shock 
mass function writes as 

N(M) = / Y[dM i N^(M i )S D (M-T[M i ), (Dl) 

where N {i \M) is the ID shock mass function along di- 
rection i. In our case all ID mass functions are iden- 
tical, iVW(M) = iV"W(M), where N (1 ~>(M) is the ID 
mass function associated with the index n studied in sec- 
tion [H21 Introducing the Mellin transform N(s) of the 
shock mass function, 



N(s) 
N(M) 



dM M s ~ N(M), 



(D2) 



c+00 ^ 



^M- s N(s), (D3) 



(where c is an arbitrary real number within the funda- 
mental strip of N(s)), we obtain at once 



N(s) = (-/V (1) (s))' 



(D4) 



Assuming that the ID mass functions show the low- 
mass power-law tails (|B3|) without logarithmic prefactors 
(which has only been proved for n = and n = —2 but 
is consistent with numerical simulations for other values 
of n, see Figs. [T] and [3]), we obtain 



n — 1 



7V (1) (s) 
N{s)~ 



1 



s + {n-l)/T 

-d 



(D5) 



, + ) . (D6) 



From standard properties of the Mellin transform this 
yields the low-mass asymptotic behavior (|75)l . Thus, we 
obtain for all dimensions d a low-mass power-law tail, 
with the same exponent (n— 1)/2, but with a logarithmic 
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prefactor to the power (d— 1). The ID high-mass cutoff 
([57]) gives the large- s behaviors (keeping only the leading- 
order term) 

s-J-oo: N^(s) ~ s s/{n+3 \ N(s) ~ s ds ^ n+3 \ (D7) 

whence the asymptotic tail ([79]) . Thus, we obtain for 
all dimensions d a high-mass modified-exponential cutoff, 
but with an exponent (n + 2>)/d that decreases at higher 
d. 

The same analysis applies to the probability distribu- 
tion, Px (t)) , of the overdensity within cubic cells of size 
X . Indeed, we again have rj = f^ i=1 f]i, where rji is the 
"ID overdensity" along direction i, and the density prob- 
ability distribution function writes as 

Px(v)= n d ^ p x (»») g o (v - n • (° 8 ) 

JO ■ 1 

i— 1 I 

This yields the Mellin transform 

Px( S ) = (pi: 1) ( S )) d . (D9) 

As for the shock mass function this leads to the high- 
density cutoff 

r/^oo: \nP x {v) - -X n+3 ^ n+3 ^ d . (D10) 

Again we recover the characteristic exponents (f74")l of the 
isotropic case. 

2. Index n = —2 
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FIG. 20: The product r] 2 Px(rj), where Px(v) is the proba- 
bility distribution function of the overdensity r\ within cubic 
cells of size X, for the separable dynamics with n = —2. Di- 
mensions d — 1 (solid line), d — 2 (dashed line), and d — 3 
(dot-dashed line), are shown for two values of the cell size X. 



For n = —2, where the ID mass function N^'(M) and 
density probability distribution PQ' (77) are given by the 



simple expressions (|B2|) and (IB5|) . we can derive explicit 
expressions. Thus, the ID Mellin transforms write as 



iV (1) (s) 



1 



(Dll) 



P 



(i 



x 



( 8 ) = 2\—e 

V 7T 



2X 



^-3/ 2 (2A),(D12) 



where K v is the modified Bessel function of the second 
kind of order v. This gives the shock mass function 
and the density probability distribution for all d through 
Eqs. (|D4|) and (|D9|) . At low dimensions it is simpler to use 
the integrals (|D1|) and (|D8|) . which gives for instance in 
2D the expressions ([5tJ ]) - (|5Tj) . This yields the asymptotic 
behaviors 



2, M -+ : N(M) 



In M 



M -3/2 



(D13) 



M -)■ 00 : N{M) - -L M" 7 / 4 e~ 2v/ ^, (D14) 



and 



which agree with the general results ([78]) . ([79]) . and 

(EES). 

We show our results in Figs. |9] and [20] for the shock 
mass function N(M) and the probability distribution 
Px(f]), in dimensions d = 1,2 and 3 and for the index 
n = — 2. To see more clearly the intermediate power- 
law regime we plot the product rj 2 Px(r)) m Fig. [201 In 
agreement with Eq. (ID10[ ). for higher d the density proba- 
bility distribution shows smoother high- and low-density 
cutoffs and a broader peak, which again can be under- 
stood as the result of a multiplicative process. As for the 
isotropic case, at smaller scales an intermediate power- 
law regime develops, as can be seen explicitly in Eq. (|8Tj) 
in 2D. However, because of logarithmic prefactors it is 
more difficult to see it in the figure for higher d. 
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